search

TinkerTools / tinker

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

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

Closed jacekkozuch closed 4 years ago

jacekkozuch commented 4 years ago

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

  2. 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 commented 4 years ago

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 commented 4 years ago

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 ???

1https://github.com/TinkerTools/tinker/issues/1 0x7f961e7fadd0 in save_parsed_format

at ../../../libgfortran/io/format.c:146

2https://github.com/TinkerTools/tinker/pull/2 0x7f961e80a187 in data_transfer_init

at ../../../libgfortran/io/transfer.c:2794

3https://github.com/TinkerTools/tinker/issues/3 0x4c83ac in ???

4https://github.com/TinkerTools/tinker/pull/4 0x404712 in ???

5https://github.com/TinkerTools/tinker/pull/5 0x414a21 in ???

6https://github.com/TinkerTools/tinker/pull/6 0x40269c in ???

7https://github.com/TinkerTools/tinker/issues/7 0x7f961d8b9494 in ???

8https://github.com/TinkerTools/tinker/pull/8 0x4026eb in ???

9https://github.com/TinkerTools/tinker/pull/9 0xffffffffffffffff in ???

Segmentation fault

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/TinkerTools/tinker/issues/67, or unsubscribehttps://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 commented 4 years ago

Hi Jay,

Thanks, I'll try that as You suggested.

Best, Jacek

On Wed, 14 Oct 2020, 22:58 Jay Ponder, notifications@github.com 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 https://github.com/TinkerTools/tinker/issues/67#issuecomment-708656133, or unsubscribe https://github.com/notifications/unsubscribe-auth/ARMBZJ34D7X4EFFHWIAHJALSKYGGHANCNFSM4SQ4YHYA .

jayponder commented 4 years ago

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.