TinkerTools / tinker

Tinker: Software Tools for Molecular Design
https://dasher.wustl.edu/tinker/
Other
130 stars 61 forks source link

NVT cannot proceed #73

Closed YangmeiLi closed 3 years ago

YangmeiLi commented 3 years ago

Dear Developers,

I am a new user to TINKER. I would like to calculate the IR spectrum of the peptide under the polarizable force field AMOEBA (amoebapro13.prm). Instead of creating the solute and inserting it into a water solvent, I created the xyz files of peptide, water box and counter ions separately and then concatenated them together. I assume the atom contacts are acceptable, as the atom coordinates are from an equilibrated system I simulated before with GROMACS under classical force fields.

The energy minimization seemed ok. But I got an energy burst in the NVT process and the program terminated. I can't really figure out why. Could you help me on this? Many thanks.

Attached are the log file, the coordinate files and the parameter files.

All the best, Yangmei tinker_files.zip

jayponder commented 3 years ago

Hi Yangmei,

There are several problems with your calculation setup in the *.key files.

(1) You are using the amoebapro13 force field. This force field uses a specific set of functional forms, and you should not change those as the parameters will be broken. For example, amoeba force fields do not use Lennard-Jones for vdw, and you can't just substitute in the L-J functional form. Similar for the polarization model, where amoeba is not intended for use with "direct" polarization. So those keywords must be removed.

(2) During the dynamics you are using very strong restraints on a set of about 16 atomic to keep them at their original positions. The default force constant for the restrain-position keyword is 100, and you have used 1000. Very stiff restraints that are still part of the dynamics are always a bad idea, and if stiff enough will make the dynamics unstable. If you change back to the default force constant of 100, then at least the dynamics will run. (Also, if you use such restraints they must be present during minimization prior to dynamics..).

A better idea, if you want those atoms to not move, is to remove their motion from the dynamics entirely. You do this in Tinker by making those atoms "inactive", which means they will be part of the energy and force calculation but won't be allowed to move during minimization or dynamics. See the suggested keyfile below.

(3) Finally, the calculations will be much more efficient if you include some additional keywords. The neighbor-list keyword will make things run much faster. With amoeba force field you should use the RESPA integrator, which will allow you to use 2.0 fs MD steps. There are also some polarization keywords that will help efficiency. You do not need to set the OpenMP threads, as by default Tinker will use all available cores on your machine. Similarly you do not need to set the archive keyword, as it is the default in recent versions of Tinker. We recommend the BUSSI thermostat over BERENDSEN, as it is just as robust, but unlike BERENDSEN it gives the correct canonical ensemble.

(4) While you are using a recent version of Tinker, we always recommend getting the latest version on Github, which will be somewhat more recent than what you are using.

Below is the *.key file that I would suggest for both minimization and dynamics. Note that I've switched to using inactive instead of restrain-position (which I've commented out, with the force constant reduced to 100).

parameters ../params/amoebapro13.prm verbose

a-axis 37.026
b-axis 36.987
c-axis 32.709

integrator RESPA thermostat BUSSI

neighbor-list vdw-cutoff 9.0

ewald ewald-cutoff 7.0

polar-eps 0.00001 polar-predict

inactive 2 18 37 49 62 78 97 109 122 inactive 138 157 169 182 198 217 229

restrain-position 2 100.0

restrain-position 18 100.0

restrain-position 37 100.0

restrain-position 49 100.0

restrain-position 62 100.0

restrain-position 78 100.0

restrain-position 97 100.0

restrain-position 109 100.0

restrain-position 122 100.0

restrain-position 138 100.0

restrain-position 157 100.0

restrain-position 169 100.0

restrain-position 182 100.0

restrain-position 198 100.0

restrain-position 217 100.0

restrain-position 229 100.0

I just ran a quick test, and your system now seems to minimize and then run dynamics correctly. So I'm going to close this issue. If you have further or new questions, please send them to my email at ponder@dasher.wustl.edu

YangmeiLi commented 3 years ago

Dear Dr. Ponder,

Thanks for your prompt reply. I really appreciate your detailed explanation and the *.key file, which are super helpful.

Cheers, Yangmei