Issue with unified AMOEBA (uAMOEBA) water model in Tinker 8.7

[jacekkozuch]

I've encountered some issues with the parameters of uAMOEBA water published in Qi et al., 2015. J. Chem. Phys., 143, 014504. Since this model is not implemented into Tinker 8.7 by default (that's the version I am currently using), I modified water03.prm using the parameters shown in Table 1 of the original paper. However, when starting minimize.x, dynamic.x or analyze.x I am observing that Tinker 8.7 gets stuck or crashes. I tried to fix these issues and isolated this issue to be due to the multipole parameters of the H atoms being zero throughout (1=O, 2=H):

multipole 2 1 2 0.00000
0.00000 0.00000 0.00000
0.00000
0.00000 0.00000
0.00000 0.00000 0.00000

Until now, I was not able to solve these issues in Tinker 8.7. I solely managed to get it to work in Tinker 6.3.3, when changing at least one multipole component to small values such as 0.000001. I attached my current uwater.prm file, that works in Tinker 6.3.3 below (with those small multipole components where water03 had some unequal to 0). Maybe, upon closer inspection, you'll find some issues that I've overlooked. Also, below the prm file are some error messages I obtained in Tinker 8.7 and 6.3.3, when all multipole components of H are zero.

Thanks for looking into this issue.

uwater.prm:
##############################
## ##
## Force Field Definition ##
## ##
##############################

forcefield uAMOEBA-WATER
#Added based on water03 and the paper cited below

bond-cubic -2.55
bond-quartic 3.793125
angle-cubic -0.014
angle-quartic 0.000056
angle-pentic -0.0000007
angle-sextic 0.000000022
opbendtype ALLINGER
opbend-cubic -0.014
opbend-quartic 0.000056
opbend-pentic -0.0000007
opbend-sextic 0.000000022
torsionunit 0.5
vdwtype BUFFERED-14-7
radiusrule CUBIC-MEAN
radiustype R-MIN
radiussize DIAMETER
epsilonrule HHG
dielectric 1.0
polarization MUTUAL
vdw-12-scale 0.0
vdw-13-scale 0.0
vdw-14-scale 1.0
vdw-15-scale 1.0
mpole-12-scale 0.0
mpole-13-scale 0.0
mpole-14-scale 0.4
mpole-15-scale 0.8
polar-12-scale 0.0
polar-13-scale 0.0
polar-14-scale 1.0
polar-15-scale 1.0
polar-12-intra 0.0
polar-13-intra 0.0
polar-14-intra 0.5
polar-15-intra 1.0
direct-11-scale 0.0
direct-12-scale 1.0
direct-13-scale 1.0
direct-14-scale 1.0
mutual-11-scale 1.0
mutual-12-scale 1.0
mutual-13-scale 1.0
mutual-14-scale 1.0

  #############################
  ##                         ##
  ##  Literature References  ##
  ##                         ##
  #############################

R. Qi, L.-P. Wang, Q. Wang, V. S. Pande, P. Ren, "United
polarizable multipole water model for molecular mechanics
simulations", J. Chem. Phys., 143, 014504 (2015)

  #############################
  ##                         ##
  ##  Atom Type Definitions  ##
  ##                         ##
  #############################

atom 1 1 O "AMOEBA Water O" 8 15.995 2
atom 2 2 H "AMOEBA Water H" 1 1.008 1

  ################################
  ##                            ##
  ##  Van der Waals Parameters  ##
  ##                            ##
  ################################

vdw 1 3.7553 0.1420
vdw 2 0.0000 0.0000 0.000

  ##################################
  ##                              ##
  ##  Bond Stretching Parameters  ##
  ##                              ##
  ##################################

bond 1 2 557.55 0.9499

  ################################
  ##                            ##
  ##  Angle Bending Parameters  ##
  ##                            ##
  ################################

angle 2 1 2 48.93 105.95

  ###############################
  ##                           ##
  ##  Urey-Bradley Parameters  ##
  ##                           ##
  ###############################

ureybrad 2 1 2 -9.74 1.5168

  ###################################
  ##                               ##
  ##  Atomic Multipole Parameters  ##
  ##                               ##
  ###################################

multipole 1 -2 -2 0.00000
0.00000 0.00000 0.70889
2.08011
0.00000 -2.08828
0.00000 0.00000 0.00817

multipole 2 1 2 0.00000
-0.00000001 0.00000 -0.00000001
-0.00000001
0.00000 -0.00000001
-0.00000001 0.00000 0.00000001

  ########################################
  ##                                    ##
  ##  Dipole Polarizability Parameters  ##
  ##                                    ##
  ########################################

polarize 1 1.7209 0.390

A few error messages that I obtained:

1. Tinker gets stuck and finally crashes, when running MINIMIZE:

SD Iter F Value G RMS F Move X Move Angle FG Call Comment
0 -8188.6155 NaN 1

2. With dynamic it starts loading parameters and ends with:
   ...
   Valence 9719-H 2 0 1
   Valence 9720-H 2 0 1
   Valence 9721-O 1 0 2
   Valence 9722-H 2 0 1
   Valence 9723-H 2 0 1
   Valence 9724-O 1 0 2
   Valence 9725-H 2 0 1
   Valence 9726-H 2 0 1
   Valence 9727-O 1 0 2
   Valence 9728-H 2 0 1
   Valence 9729-H 2 0 1

Tinker is Unable to Continue; Terminating the Current Calculation

3. Though it's possible to start analyze, using option M I get a segmentation fault:

Total Electric Charge : 0.00000 Electrons

Dipole Moment Magnitude : 127.490 Debye

Dipole X,Y,Z-Components : -34.875 105.286 62.868

Quadrupole Moment Tensor : -246.924 -879.025 1794.891
(Buckinghams) -879.025 -160.170 -1411.151
1794.891 -1411.151 407.095

Principal Axes Quadrupole : -1801.976 -1016.433 2818.409
Radius of Gyration : 22.992 Angstroms

Center of Mass Coordinates : -0.111634 0.048246 0.053490
Euler Angles (Phi/Theta/Psi) : -64.611 -22.599 149.178

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0 0x7f961d8cd33f in ???
#1 0x7f961e7fadd0 in save_parsed_format
at ../../../libgfortran/io/format.c:146
#2 0x7f961e80a187 in data_transfer_init
at ../../../libgfortran/io/transfer.c:2794
#3 0x4c83ac in ???
#4 0x404712 in ???
#5 0x414a21 in ???
#6 0x40269c in ???
#7 0x7f961d8b9494 in ???
#8 0x4026eb in ???
#9 0xffffffffffffffff in ???
Segmentation fault

[jayponder]

Thanks for sending this. A "workaround" that avoids the problem is to add a tiny non-zero polarizability parameter for the water hydrogen. I suggest adding a line like this to your parameter file:

polarize 2 0.00000001 0.39

Note the tiny value is the atomic polarizability of the hydrogen. The other number, 0.39, is Thole damping value used with AMOEBA-like models, and should be left at 0.39. Also, I would note you should be able to set all the hydrogen multipole components to exactly zero.

Obviously, this is a "bug", and is almost certainly due to the fact that Tinker tries to remove atoms with no polarizability (either exactly zero or not defined) from the polarization calculation. This requires an indexing scheme between the full set of atoms and the sites with non-zero polarizability. This indexing is likely broken. But if you just make sure all atoms have at least some polarizability then things work correctly.

I will look into this further, and fix the bug. But in the meantime, the suggested workaround seems safe.

[pren]

Jacek, I also think you should use all-atom AMOEBA water instead unless you have compelling reason to use uamoeba. The GPU code should be plenty fast. Thanks Pengyu From: Jacek Kozuch <notifications@github.com> Sent: Wednesday, October 14, 2020 12:25 PM To: TinkerTools/tinker <tinker@noreply.github.com> Cc: Subscribed <subscribed@noreply.github.com> Subject: [TinkerTools/tinker] Issue with unified AMOEBA (uAMOEBA) water model in Tinker 8.7 (#67) I've encountered some issues with the parameters of uAMOEBA water published in Qi et al., 2015. J. Chem. Phys., 143, 014504. Since this model is not implemented into Tinker 8.7 by default (that's the version I am currently using), I modified water03.prm using the parameters shown in Table 1 of the original paper. However, when starting minimize.x, dynamic.x or analyze.x I am observing that Tinker 8.7 gets stuck or crashes. I tried to fix these issues and isolated this issue to be due to the multipole parameters of the H atoms being zero throughout (1=O, 2=H): multipole 2 1 2 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 Until now, I was not able to solve these issues in Tinker 8.7. I solely managed to get it to work in Tinker 6.3.3, when changing at least one multipole component to small values such as 0.000001. I attached my current uwater.prm file, that works in Tinker 6.3.3 below (with those small multipole components where water03 had some unequal to 0). Maybe, upon closer inspection, you'll find some issues that I've overlooked. Also, below the prm file are some error messages I obtained in Tinker 8.7 and 6.3.3, when all multipole components of H are zero. Thanks for looking into this issue. uwater.prm: ############################## ## ## ## Force Field Definition ## ## ## ############################## forcefield uAMOEBA-WATER #Added based on water03 and the paper cited below bond-cubic -2.55 bond-quartic 3.793125 angle-cubic -0.014 angle-quartic 0.000056 angle-pentic -0.0000007 angle-sextic 0.000000022 opbendtype ALLINGER opbend-cubic -0.014 opbend-quartic 0.000056 opbend-pentic -0.0000007 opbend-sextic 0.000000022 torsionunit 0.5 vdwtype BUFFERED-14-7 radiusrule CUBIC-MEAN radiustype R-MIN radiussize DIAMETER epsilonrule HHG dielectric 1.0 polarization MUTUAL vdw-12-scale 0.0 vdw-13-scale 0.0 vdw-14-scale 1.0 vdw-15-scale 1.0 mpole-12-scale 0.0 mpole-13-scale 0.0 mpole-14-scale 0.4 mpole-15-scale 0.8 polar-12-scale 0.0 polar-13-scale 0.0 polar-14-scale 1.0 polar-15-scale 1.0 polar-12-intra 0.0 polar-13-intra 0.0 polar-14-intra 0.5 polar-15-intra 1.0 direct-11-scale 0.0 direct-12-scale 1.0 direct-13-scale 1.0 direct-14-scale 1.0 mutual-11-scale 1.0 mutual-12-scale 1.0 mutual-13-scale 1.0 mutual-14-scale 1.0 ############################# ## ## ## Literature References ## ## ## ############################# R. Qi, L.-P. Wang, Q. Wang, V. S. Pande, P. Ren, "United polarizable multipole water model for molecular mechanics simulations", J. Chem. Phys., 143, 014504 (2015) ############################# ## ## ## Atom Type Definitions ## ## ## ############################# atom 1 1 O "AMOEBA Water O" 8 15.995 2 atom 2 2 H "AMOEBA Water H" 1 1.008 1 ################################ ## ## ## Van der Waals Parameters ## ## ## ################################ vdw 1 3.7553 0.1420 vdw 2 0.0000 0.0000 0.000 ################################## ## ## ## Bond Stretching Parameters ## ## ## ################################## bond 1 2 557.55 0.9499 ################################ ## ## ## Angle Bending Parameters ## ## ## ################################ angle 2 1 2 48.93 105.95 ############################### ## ## ## Urey-Bradley Parameters ## ## ## ############################### ureybrad 2 1 2 -9.74 1.5168 ################################### ## ## ## Atomic Multipole Parameters ## ## ## ################################### multipole 1 -2 -2 0.00000 0.00000 0.00000 0.70889 2.08011 0.00000 -2.08828 0.00000 0.00000 0.00817 multipole 2 1 2 0.00000 -0.00000001 0.00000 -0.00000001 -0.00000001 0.00000 -0.00000001 -0.00000001 0.00000 0.00000001 ######################################## ## ## ## Dipole Polarizability Parameters ## ## ## ######################################## polarize 1 1.7209 0.390 A few error messages that I obtained: 1. Tinker gets stuck and finally crashes, when running MINIMIZE: SD Iter F Value G RMS F Move X Move Angle FG Call Comment 0 -8188.6155 NaN 1 1. With dynamic it starts loading parameters and ends with: ... Valence 9719-H 2 0 1 Valence 9720-H 2 0 1 Valence 9721-O 1 0 2 Valence 9722-H 2 0 1 Valence 9723-H 2 0 1 Valence 9724-O 1 0 2 Valence 9725-H 2 0 1 Valence 9726-H 2 0 1 Valence 9727-O 1 0 2 Valence 9728-H 2 0 1 Valence 9729-H 2 0 1 Tinker is Unable to Continue; Terminating the Current Calculation 1. Though it's possible to start analyze, using option M I get a segmentation fault: Total Electric Charge : 0.00000 Electrons Dipole Moment Magnitude : 127.490 Debye Dipole X,Y,Z-Components : -34.875 105.286 62.868 Quadrupole Moment Tensor : -246.924 -879.025 1794.891 (Buckinghams) -879.025 -160.170 -1411.151 1794.891 -1411.151 407.095 Principal Axes Quadrupole : -1801.976 -1016.433 2818.409 Radius of Gyration : 22.992 Angstroms Center of Mass Coordinates : -0.111634 0.048246 0.053490 Euler Angles (Phi/Theta/Psi) : -64.611 -22.599 149.178 Program received signal SIGSEGV: Segmentation fault - invalid memory reference. Backtrace for this error: #0 0x7f961d8cd33f in ??? #1<#1> 0x7f961e7fadd0 in save_parsed_format at ../../../libgfortran/io/format.c:146 #2<#2> 0x7f961e80a187 in data_transfer_init at ../../../libgfortran/io/transfer.c:2794 #3<#3> 0x4c83ac in ??? #4<#4> 0x404712 in ??? #5<#5> 0x414a21 in ??? #6<#6> 0x40269c in ??? #7<#7> 0x7f961d8b9494 in ??? #8<#8> 0x4026eb in ??? #9<#9> 0xffffffffffffffff in ??? Segmentation fault — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#67>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABNC6XRZLB3XFFGJNY4QLJDSKXNHZANCNFSM4SQ4YHYA>. This message is from an external sender. Learn more about why this matters.<https://ut.service-now.com/sp?id=kb_article&number=KB0011401>

[jacekkozuch]

Hi Jay, Thanks, I'll try that as You suggested. Best, Jacek

…

On Wed, 14 Oct 2020, 22:58 Jay Ponder, ***@***.***> wrote: Thanks for sending this. A "workaround" that avoids the problem is to add a tiny non-zero polarizability parameter for the water hydrogen. I suggest adding a line like this to your parameter file: polarize 2 0.00000001 0.39 Note the tiny value is the atomic polarizability of the hydrogen. The other number, 0.39, is Thole damping value used with AMOEBA-like models, and should be left at 0.39. Also, I would note you should be able to set all the hydrogen multipole components to exactly zero. Obviously, this is a "bug", and is almost certainly due to the fact that Tinker tries to remove atoms with no polarizability (either exactly zero or not defined) from the polarization calculation. This requires an indexing scheme between the full set of atoms and the sites with non-zero polarizability. This indexing is likely broken. But if you just make sure all atoms have at least some polarizability then things work correctly. I will look into this further, and fix the bug. But in the meantime, the suggested workaround seems safe. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#67 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ARMBZJ34D7X4EFFHWIAHJALSKYGGHANCNFSM4SQ4YHYA> .

[jayponder]

I believe this is now fixed in the sources here on Github, and the workaround above is no longer needed if using the updated Github sources. It will also be fixed on https://dasher.wustl.edu/tinker/ the next time the version on that site is updated. I've also added a new uwater.prm parameter file in the /params area of the Tinker distribution that implements this water model.