Question about protein-ligand simulation

[xtzhang0216]

lig.key.txt

Hello Tinker developers,

I am trying to set up a protein–ligand molecular dynamics simulation in Tinker using the AMOEBA polarizable force field. I recently transitioned from GROMACS to Tinker specifically for AMOEBA, so I would like to confirm that my workflow and parameter handling are correct.

Background

• Protein force field: amoebabio18.prm

• Ligand parameters: generated using poltype2 (lig.xyz and lig.key)

• Initial structure: a protein–ligand–counterion complex that was previously equilibrated using GROMACS with the CHARMM force field, saved as a PDB file

• My goal is to convert this structure into Tinker format and then run AMOEBA MD (minimization → equilibration → production)

---

Issue 1: Loading protein and ligand parameters

After generating the ligand parameters with poltype2, I attempted to convert the PDB structure into a Tinker XYZ file using pdbxyz.

I created a complex.key file with the following contents:

parameters amoebabio18.prm
parameters lig.key

Then I ran:

pdbxyz.x complex.pdb -k complex.key

However, the output XYZ file is empty, which is confusing to me. I would appreciate clarification on the correct and supported way to load ligand parameters together with amoebabio18.prm.

---

Issue 2: Recommended workflow for solvation and MD with AMOEBA

After the protein–ligand system is correctly loaded, I would also like to confirm the recommended workflow for AMOEBA simulations in Tinker.

My current understanding is roughly:

1. Solvation

   • Use AMOEBA water (or TIP4P or else? Please let me know if any other water model is recommended), to create a water box.

   • Fit the complex coordinate into the water box.

2. Simulation

   • Energy minimization

   • Equilibration

   • Production MD

Could you please confirm whether this workflow is correct?

Additionally, if possible, I would greatly appreciate:

• A template example key file for minimization/equilibration/production with amoeba forcefield

• Especially one that includes heavy-atom positional restraints for equilibration and production

---

I apologize for a lot of questions, but I am still learning the Tinker/AMOEBA workflow and want to make sure I am using it correctly.
Thank you very much for your time and for maintaining Tinker—any advice or clarification would be extremely helpful.

Best regards,
Xintong

---

[jayponder]

Unfortunately, some of this process is perhaps not as easy as it should be for first time users of Tinker. And it's rather different from the way CHARMM, Amber, GROMACS, etc. operate.

Issue 1:

You should convert the protein and your ligand to Tinker XYZ format separately, then combine them into one file. So make a PDB file with the protein, and convert this file to Tinker XYZ format. For example, if you have a file named "protein.pdb", then make "protein.key" with the line "parameters amoebabio18.prm". (Note that if the amoebabio18.prm file is not in the directory you are working in, then the parameters keyword in the .key file needs the full path name to the .prm file. Run the "pdbyz" program giving "protein.pdb" as the input file, which should create protein.xyz and protein.seq files with the coordinates in Tinker format and the sequence, respectively. You can then check the output file by running the "analyze" program using the "protein.xyz" file as input. This should report an energy for the protein without any warnings or errors.

You have already used Poltype to generate "lig.xyz" and "lig.key" files for your ligand. Again, you can use the "analyze" program to check that the "lig.xyz" file is correct and valid. Note it is also important that the atom type numbers for the ligand atoms do not conflict with the atom type numbers used for the protein in the amoebabio18 parameter file. Poltype should have used higher atom type numbers for the ligand by default.

You can then merge the "protein.xyz" file and "lig.xyz" files into a single Tinker coordinates file. Us the "xyzedit" program and run it with "protein.xyz" as input. Choose the xyzedit option to "Append a Second XYZ File to the Current One", and give the "lig.xyz" file as the second xyz file. This should produce a new version of the protein file (ie, "protein.xyz_2") that has the ligand appended below the protein coordinates. You can rename the "protein.xyz_2" to be, for example, "complex.xyz". Then also copy "lig.key" to be "complex.key", as this .key file should point to the amoebabio18 parameter file and also contain the ligand parameters for the ligand within the .key file. Once again, you can use the "analyze" program with "complex.xyz" as input to check for correctness.

Note that only one parameter (.prm) file can be included in the .key file for a calculation. Any other paraemters you need, beyond the amoebabio18 values, must be included directly in the .key file. This is something that we should change in a future release of Tinker- ie, to allow the .key file to point to multiple parameter files. So the "complex.key" file you describe above with two parameters lines, one for amoebabio18.prm and one for lig.key will not work.

Issue 2:

Once you have a valid Tinker files for the combined protein and ligand, for example "complex.xyz" and "complex.key", you can solvate this complex. First, you should translate the complex coordinates so the system is centered at (0,0,0). This can be done by running the "xyzedit" program with "complex.xyz" as input and using the option to "Translate Center of Mass to the Origin". This will create a new output file as "complex.xyz_2".

This centered complex can then be placed into a water box. You must use AMOEBA water if you are using AMOEBA parameters for the protein-ligand complex. We provide several different sizes of equilibrated AMOEBA water boxes on the Ponder lab website at https://dasher.wustl.edu/tinker/. Ideally, pick a water box with sides 20-30 Angstroms larger than the maximum interatomic distance in your protein-ligand complex. First you will need to change the atom type numbers in the equilibrated water box, as they are supplied with water oxygen as atom type 1 and hydrogen at atom type 2. In the amoebabio18 parameters, AMOEBA water oxygen is atom type 349 and the hydrogens are atom type 350. You can manually change these atom type numbers in the water box .xyz file, or they can be changed by using the "xyzedit" program and the option to "Replace Old Atom Type with a New Type". Then run the "xyzedit" program with "complex.xyz" as the input file, and use the option to "Soak Current Molecule in a Box of Solvent", giving the name of the AMOEBA water box with modified atom types when prompted for the solvent. This will place the complex in the water box, and remove water that overlaps the complex. If you wish to add ions, this can be done using yet another "xyzedit" option to "Place Monoatomic Ions around a Solute", which will remove water molecules outside a buffer zone around your complex and replace them with your specified number and type of ions from the amoebabio18 parameters.

As you suggest, you will then want to minimize the final solvated complex, run some equilibration MD and then production MD. First you will want to add keywords to your "complex.key" to turn on various options for efficient minimization and molecular dynamics- things like PME, use of neighbor lists, etc. Several sample .key file setups are provided in the /bench and /example areas of the Tinker distribution. For example, you might look at the bench7.key file in the /bench benchmark directory.

For minimization, you should use the Tinker "minimize" program. You need only a short minimization to loose convergence to prepare for MD. I would suggest using an RMS gradient cutoff of 1.0 or 3.0 instead of the default value of 0.01, which will be very slow and unnecessary. You can then use the "dynamic" program to run MD, following whatever protocols you have used with other force fields for length of equilibration, ramping temperature, etc. Tinker and AMOEBA behave rather like other packages and force fields in this regard. Since AMOEBA is a fully flexible force field, we do recommend using a RESPA multiple time step integrator (keyword: integrator RESPA) as this allows use of 2 fs time steps with AMOEBA (or even 3 fs if you turn on hydrogen mass redistribution via keyword: heavy-hydrogen). Standard, simple velocity Verlet or Beeman integrators must be run using 1.0 fs steps, or better with 0.5 fs steps, when using AMOEBA.

Finally I would note that the CPU version of Tinker is very, very slow for production MD! So for actual production simulations, you will almost certainly want to switch to the GPU-enabled Tinker-GPU or Tinker-HP codes to run on NVIDIA GPUs. These GPU codes are file compatible with the CPU Tinker code, so you can prepare files and test all of the above with the CPU code, then switch to the much faster GPU code for production simulation.

[xtzhang0216]

Thank you very much for your prompt and patient reply. I feel that I am now very close to having a fully working setup, but I am still encountering a few issues that I would appreciate your advice on.

System preparation and box shape

My protein is roughly rectangular in shape, and during equilibration I plan to apply positional restraints on heavy atoms. To reduce the number of solvent molecules and improve efficiency, I would therefore prefer to use a rectangular periodic box rather than a cubic one. I noticed that the pre-equilibrated AMOEBA water boxes provided with Tinker appear to be cubic, which makes them difficult to reuse directly in this case. As a workaround, I used a structure that had already been equilibrated under GROMACS + CHARMM, containing protein, water, and counterions, and converted it to Tinker format. I was able to:

• successfully convert the system to Tinker XYZ format,

• merge it with the POLTYPE-generated ligand XYZ,

• and obtain a combined complex.xyz and complex.key.

At this stage, the system appears to be correctly assembled.

---

Minimization: CPU vs GPU behavior

When I perform energy minimization using the CPU version, the behavior appears reasonable, with the energy decreasing monotonically:

minimize.x complex.xyz -k min.key

Enter RMS Gradient per Atom Criterion [0.01] : 1

Limited Memory BFGS Quasi-Newton Optimization :
QN Iter  F Value     G RMS     F Move     X Move   Angle  FG Call  Comment
0   0.1685D+11 0.8013D+07
1   0.1462D+11 0.7114D+07 0.2231D+10 0.0053   0.00     4  Success
2   0.1206D+11 0.6234D+07 0.2560D+10 0.0071  11.52     8  Success
3   0.1033D+11 0.5551D+07 0.1729D+10 0.0053   1.70    11  Success
4   0.8797D+10 0.4919D+07 0.1536D+10 0.0053   1.53    14  Success

However, when I use the GPU version of minimization:

minimize9 complex.xyz -k minim.key

I observe a different behavior:

Limited Memory BFGS Quasi-Newton Optimization :
QN Iter  F Value     G RMS
0  -0.2141D+10 0.2239D+08

LBFGS  --  Normal Termination due to SmallFct

Final Function Value :  -2140959388.61276102
Final RMS Gradient   :   22387290.71714260
Final Gradient Norm  : 5260012940.21560860

The GPU minimization terminates very quickly, with a different total energy and a relatively large RMS gradient. I am unsure whether this behavior is expected or indicates that the GPU minimization did not adequately relax the system.

Because the CPU minimization is very slow for this system, I proceeded using the GPU-minimized structure as the starting point for equilibration.

---

NPT equilibration issues (GPU vs CPU)

When I attempt to run NPT equilibration using the GPU version:

dynamic9 min.xyz -k NPT.key

the simulation terminates with an induced dipole convergence error:

 Enter the Number of Dynamics Steps to be Taken :  10000

 Enter the Time Step Length in Femtoseconds [1.0] :  

 Enter Time between Saves in Picoseconds [0.1] :  

 Available Statistical Mechanical Ensembles :
    (1) Microcanonical (NVE)
    (2) Canonical (NVT)
    (3) Isoenthalpic-Isobaric (NPH)
    (4) Isothermal-Isobaric (NPT)
 Enter the Number of the Desired Choice  [1] :  4

 Enter the Desired Temperature in Degrees K [298] :  

 Enter the Desired Pressure in Atm [1.0] :  

 GPU Device :  Setting Device ID to 0 from GPU utilization

 Bussi Thermostat
 Temperature               298.0 Kelvin
 Tau-Temperature           0.200 ps

 Monte Carlo Barostat
 Volume Move               100.0 Ang**3
 Pressure                  1.000 Atm
 Tau-Pressure              2.000 ps
 Volume Trial                 25
 Isotropic Fluctuation

 NRespa                        4

 Molecular Dynamics Trajectory via r-RESPA MTS Algorithm
 Backtrace
    1  dynamic9                                                      0x565b39
    2  dynamic9                                                      0x430192
    3  dynamic9                                                      0x559081
    4  dynamic9                                                      0x41b0c3
    5  dynamic9                                                      0x5948a9
    6  dynamic9                                                      0x562e92
    7  dynamic9                                                      0x563547
    8  dynamic9                                                      0x5a48ca
    9  dynamic9                                                      0x56fca7
   10  dynamic9                                                      0x58b74d
   11  dynamic9                                                      0x40d335
   12  /lib64/libc.so.6                                              __libc_start_main
   13  dynamic9                                                      0x40debe
 Terminating with uncaught exception :  INDUCE  --  Warning, Induced Dipoles are not Converged at /tmp/tinker-gpu/src/cu/amoeba/pcg.cu:183

[input.zip](https://github.com/user-attachments/files/24292743/input.zip)

In contrast, when I attempt NPT equilibration using the CPU version, the run fails with:

GRADIENT  --  Illegal Value for the Total Potential Energy

This makes me wonder whether:

1. A potential problem with my CPU or GPU build of Tinker. Both versions were built from the latest GitHub commit using the standard build procedure, without any local code modifications.

2. The initial structure converted from CHARMM (protein, water, and ions) may not be fully compatible with AMOEBA without additional preparation.

3. My simulation protocol and/or key file settings (for minimization or early NPT equilibration) may need to be modified.

---

To facilitate debugging, I have bundled and uploaded all relevant input files. Thank you again for your time and help. Any advice would be extremely valuable as I continue learning the recommended AMOEBA workflow in Tinker.

Best regards,
Xintong

[xtzhang0216]

input.zip

[jayponder]

The protein portion of your input file looks correct. However, the ligand structure is badly damaged. The ligand has two bonded atoms that are 52 Angstroms apart, many atoms that are much too close together, and other geometry problems. It does look like most of the ligand atoms are located in a pocket near the surface of the protein, so at least the general location of the ligand within the system may be correct, I'm not sure. I can tell this by viewing the complex.xyz file you sent in the Tinker-FFE GUI for Tinker. If you don't have that GUI, then you should be able to convert the complex.xyz file back into a PDB file (you will need the .seq file created when you originally converted the protein from PDB to Tinker xyz, and have it named complex.seq). You could then view the complex.pdb in VMD or any other program that views PDB files. You can also run the Tinker "analyze" program and it will output the very large interactions. For example, if you use "analyze" it shows the extremely long bond is between atoms 55173 and 55185.

So, first you need to get the ligand structure correct and placed in the binding pocket in the protein to use as your complex.xyz file. You cannot expect minimization to fix a crazy structure with 50 Angstrom bonds, etc.

Also, in the complex.key file, you should turn on particle mesh Ewald (PME) and use of neighbor lists before doing any Tinker calculations on the system, such as "analyze", running minimization or doing MD. Note that unlike some modeling codes, these things are not turned on by default in Tinker. Add the following lines to your complex.key file:

EWALD
EWALD-CUTOFF 7.0
NEIGHBOR-LIST

Including those lines will make even the CPU calculations MUCH faster!

As you may know, it can be dangerous to use a "rectangular" box and restraints just to reduce the number of atoms in you system. It is your choice, but I would probably use a larger cubic water box without restraints, perhaps 95 or 100 Angstroms. Alternatively Tinker supports both truncated octahedon and rhombic dodecahedron periodic box shapes. Just a suggestion!

Also, it looks like your protein is "capped" at the N-terminal and C-terminal ends with ACE and NME capping residues, respectively. Of course, native proteins generally do not have such caps. Are you sure this is what you want for the termini of your protein?

[xtzhang0216]

Sorry for my careless mistake earlier. I have now corrected the ligand coordinates, and I used analyze.x with the “List of the Large Individual Interactions” option to confirm that there are no obvious topology problems. After that, the energy minimization finished smoothly. Thank you very much for your patient and generous help.

I still have one question regarding the simulation setup.
When I apply position restraints on heavy atoms (force constant = 2.5), I encounter the error:

Warning, Induced Dipoles are not Converged

However, if I do not apply the restraints, the simulation runs normally.
In this case, is setting those atoms as inactive the only practical solution? I would prefer not to do so, since I still want to keep some flexibility in the system.

inputs.zip

In addition, I would like to ask a more general question. When using the AMOEBA force field, is there a recommended way to analyze protein–ligand binding free energy? For example, is it acceptable to convert the system to an AMBER or GROMACS topology for approximate analyses such as MM/PBSA or MM/GBSA? I understand that such conversions may not be fully consistent with a polarizable force field, so any guidance or best practices would be greatly appreciated.

Thank you again for your help.

[xtzhang0216]

Hello developer, and thank you very much for your help so far.

I have one additional question regarding performance. I am running the latest tinker-gpu on an RTX 4090, with a system of ~55,500 atoms and a 2 fs timestep. The observed speed is about 24 ns/day.

I would like to ask whether this performance is reasonable for AMOEBA, and whether other engines such as OpenMM are expected to be faster for systems of this size, or if this is already close to optimal on consumer GPUs.

Thank you very much for your time, and Merry Christmas 🎄!

[jayponder]

To compute ligand binding free energies, we have used a fairly standard windowed method in lambda, and then the BAR procedure to extract the free energies. This is documented in several papers over the last 15+ years from my group at Washington University and that of Pengyu Ren at Univ Texas-Austin. The "bar" program that is part of Tinker and Tinker-GPU can be used to extract the free energies. There is no executable program we distribute to automate running the MD windows at various lambda values. But this involves just some simple changes in .key files to define the "ligand" and to set the lambda values for the vdw and electrostatic components. It is simple to write some scripts of your own to mostly automate the necessary MD windows. As an example, see "Accurate Host-Guest Binding Free Energies Using the AMOEBA Polarizable Force Field", M. K. J. Chung, R. J. Miller, B. Novak, Z. Wang and J. W. Ponder, J. Chem. Inf. Model., 63, 2769-2782 (2023), and many earlier papers.

If you want to implement an "MM/GBSA" method using AMOEBA, then you should use the corresponding GB-like implicit solvent model for AMOEBA, which is called Generalized Kirkwood (so GK, instead of GB). The GK model is turned on in the .key file via the keyword "SOLVATE GK". The GK model extends GB to handle the polarization and atomic multipoles used by AMOEBA.

Finally, I have some comments regarding speed of calculation. There are some other keywords you can use to further speed the calculation in ns/day. You can use "POLAR-PREDICT" and "POLAR-EPS 0.00001". Neither of those keywords will have any effect on accuracy of the calculation. You can also try "POLARIZATION OPT3" which changes the induced dipole model slightly, but is generally safe. Similar to simple fixed charge force fields, you can also run the RESPA integrator with AMOEBA at 3 fs time steps by doing hydrogen mass redistribution (add keyword "HEAVY-HYDROGEN"). This maintains the canonical ensemble, but does change transport properties somewhat. However it is generally safe for the sampling to compute free energies.

While I don't have an RTX 4090 to run benchmarks, I would say that the above ideas may give you somewhat more than the 24 ns/day you report for 55.5K atoms. However, I don't think you will do a large amount better even with these tricks. The OpenMM program will not give better speed for AMOEBA than Tinker-GPU, in fact I expect it will be similar or slightly slower.

[xtzhang0216]

Hi developer,

Thank you very much for your detailed reply—I really appreciate you taking the time to answer, especially during the holidays.

I just wanted to briefly follow up on one earlier question, in case I was not clear before: do you have any general recommendations regarding the use of restraints to avoid " Induced Dipoles are not Converged"?

No worries at all if you don’t have time to comment on this now. Thanks again for your help.

Best regards,
Xintong

[jayponder]

The "Induced Dipoles are not Converged" message is often the first point of failure with AMOEBA, and can be caused by many things. The most common reason is simply a very bad initial structure. Polarizable force fields in general are more sensitive to unreasonable structures than simple fixed charge force fields. Lack of convergence for the induced dipoles can also be caused by bad or extreme multipole parameters. But since you are using amoebabio18 values for the protein and water, and Poltype values for the ligand, the parameters should be good.

If you are doing an initial minimization, perhaps before starting MD, and getting the induced dipole convergence error, then there are three levels of things to try to avoid the error. First would be to use direct polarization instead of full mutual polarization by adding the "POLARIZATION DIRECT" keyword. Second would be to turn off polarization entirely via the "POLARIZETERM NONE" keyword. Finally would be to turn off atomic multipole electrostatics and polarization by using "MULTIPOLETERM NONE" in addition to "POLARIZETERM NONE". Any of the above should just be used for a first minimization with a loose convergence criterion. Then remove the above keywords and repeat the minimization with the full AMOEBA model.