TinkerTools / tinker

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

Optimize issue with TIP5P? #154

Closed mircodaniel closed 3 months ago

mircodaniel commented 4 months ago

When I try to use "optimize" tool with following .xyz file I get a strange result: HAc_H2O.xyz: 13 HAc_H2O 1 C 0.558027 0.163473 0.126991 209 2 3 5 2 O 1.232527 0.683763 1.007634 210 1 3 O 1.081084 -0.259078 -1.051402 211 1 4 4 H 2.008161 -0.080648 -1.050340 212 3 5 C -0.939677 -0.081389 0.223051 80 1 6 7 8 6 H -1.462501 0.448919 -0.572582 85 5 7 H -1.156392 -1.145769 0.135116 85 5 8 H -1.321230 0.270729 1.181533 85 5 9 O 3.680000 -0.002207 -0.072659 73 10 11 12 13 10 H 3.680000 -0.741019 0.535938 74 9 11 H 3.679998 0.772184 0.489969 74 9 12 Lp 3.108456 -0.014479 -0.476624 75 9 13 Lp 4.251544 -0.014479 -0.476624 75 9

HAc_H2O.key: PARAMETERS "tinker/params/oplsaa.prm" DIELECTRIC 1.0

HAc_H2O_afterOpt.xyz: 13 HAc_H2O 1 C 0.729027 0.194135 0.376030 209 2 3 5 2 O 1.524620 0.860076 1.328523 210 1 3 OH 1.342037 -0.449580 -0.696372 211 1 4 4 HO 2.582875 -0.059567 -0.506366 212 3 5 CT -0.845532 -0.184949 -0.272224 80 1 6 7 8 6 HC -1.460677 0.488812 -0.465849 85 5 7 HC -1.121235 -1.032908 0.053074 85 5 8 HC -1.267620 0.104566 0.680040 85 5 9 OW 3.298847 0.173735 -0.319447 73 10 11 12 13 10 HW 3.465311 -0.410305 0.365589 74 9 11 HW 3.263412 0.511460 0.393107 74 9 12 LP 2.582875 -0.059567 -0.506366 75 9 13 LP 4.306055 -0.135908 -0.429738 75 9

I assume "optimize" can't handle TIP5P water model. Is my assumption correct?

jayponder commented 3 months ago

Yes, your assumption is correct. TIP5P is intended (by Bill Jorgensen, the originator of TIP5P) to be a rigid water model. In the OPLS-AA parameter file distributed with Tinker, we added some fairly stiff bond and angle parameters so that (most) calculations trying to use "flexible TIP5P" will at least run without failure. But in your case, since the TIP5P has no vdw directly on the lone pair sites, and the carboxylic acid hydrogen also is not a vdw site, the electrostatic forces overcome the TIP5P bond and angle values we have added to OPLS-AA and you get fusion of those sites. You could "fix" this by doing minimization between rigid acetic acid and TiP5P molecules, or perhaps by making the TIP5P bond and angle parameter force constants even larger still.

Note the Tinker OPLS-AA parameter file contains several Jorgensen-developed water models: TIP3P, TIP4P, TIP5P, TIP3F and TIP4F. The intent is that TIP3P, TIP4P and TIP5P are rigid models, while TIP3F and TIP4F were calibrated with "officially approved" parameters for use as flexible models. The F versions are obscure and have never been much used. The most common application, of course, is to use the rigid versions in MD simulation via SHAKE or RATTLE to keep the water molecules rigid. To my knowledge, there is no TIP5F version. And there are many other "non-Jorgensen" derivative models in the literature- for example TIP4P-EW was reparameterized for use with Ewald summation, and has received quite a lot of use.