qTIP4P/f force field model
93 views
Skip to first unread message
Ivan Jackson
Mar 9, 2020, 2:16:39 PM3/9/20
to cp2k
Hi,
I am trying to run bulk water NVT simulation using Scott Habershon's q-TIP4P/F model (J. Chem. Phys. 131, 024501 (2009)) in CP2K. Starting from TIP4P-2005, I changed the TIP4P-2005 intermolecular parameters to be as defined for q-TIP4p/f model. However, when I run the simulation, it crashes after a few steps. When I analyse the resultant trajectory, I noticed that the OH bond length increases to 1.8 angstroms (far away from r_0=0.94 angstrom for q-TIP4P/f), where as I able to run the original TIP4P-2005 without any issues. Do you any idea about how to resolve this particular issue for the q-TIP4P/f model? If any there are any inconsistencies in my input, I would be really grateful if you could point them out. The relevant part of force field model in my input file is:
&FORCE_EVAL
METHOD FIST
&MM
&FORCEFIELD
&SPLINE
RCUT_NB 0.9000000000000000E+01
EMAX_SPLINE 1.0000000000000000E+06
&END SPLINE
&NONBONDED
&LENNARD-JONES
ATOMS OW OW
EPSILON 9.3200000000000003E+01
SIGMA 3.1589000000000000E+00
&END LENNARD-JONES
&LENNARD-JONES
ATOMS OW MW
EPSILON 0.0000000000000000E+00
SIGMA 1.0000000000000001E-01
&END LENNARD-JONES
&LENNARD-JONES
ATOMS OW HW
EPSILON 0.0000000000000000E+00
SIGMA 1.0000000000000001E-01
&END LENNARD-JONES
&LENNARD-JONES
ATOMS HW HW
EPSILON 0.0000000000000000E+00
SIGMA 1.0000000000000001E-01
&END LENNARD-JONES
&LENNARD-JONES
ATOMS HW MW
EPSILON 0.0000000000000000E+00
SIGMA 1.0000000000000001E-01
&END LENNARD-JONES
&LENNARD-JONES
ATOMS MW MW
EPSILON 0.0000000000000000E+00
SIGMA 1.0000000000000001E-01
&END LENNARD-JONES
&END NONBONDED
&CHARGE
ATOM OW
CHARGE 0.0000000000000000E+00
&END CHARGE
&CHARGE
ATOM HW
CHARGE 5.5640000000000001E-01
&END CHARGE
&CHARGE
ATOM MW
CHARGE -1.1128000000000000E+00
&END CHARGE
&BOND
ATOMS OW HW
KIND G87
K 5.3550462337743521E-01
R0 1.8088458543981376E+00
&END BOND
&BOND
ATOMS OW MW
KIND G87
K 0.0000000000000000E+00
R0 3.3070207325498757E-01
&END BOND
&BOND
ATOMS HW HW
KIND G87
K 0.0000000000000000E+00
R0 3.3070207325498757E-01
&END BOND
&BOND
ATOMS HW MW
KIND G87
K 0.0000000000000000E+00
R0 3.3070207325498757E-01
&END BOND
&BEND
ATOMS HW OW HW
KIND G87
K 0.0000000000000000E+00
THETA0 1.9111355309337910E+00
&END BEND
&BEND
ATOMS HW OW MW
KIND G87
K 0.0000000000000000E+00
THETA0 9.5556776546689548E-01
&END BEND
&BEND
ATOMS MW HW OW
KIND G87
K 0.0000000000000000E+00
THETA0 1.9111355309337910E+00
&END BEND
&BEND
ATOMS MW HW HW
KIND G87
K 0.0000000000000000E+00
THETA0 1.9111355309337910E+00
&END BEND
&BEND
ATOMS OW HW HW
KIND G87
K 0.0000000000000000E+00
THETA0 1.9111355309337910E+00
&END BEND
&BEND
ATOMS OW MW HW
KIND G87
K 0.0000000000000000E+00
THETA0 1.9111355309337910E+00
&END BEND
&BEND
ATOMS HW MW HW
KIND G87
K 0.0000000000000000E+00
THETA0 1.9111355309337910E+00
&END BEND
&END FORCEFIELD
My input file is attached. Many thanks for your help.
Best Regards.
Ivan Jackson
Thomas Kühne
Mar 9, 2020, 3:05:33 PM3/9/20
to 'Dorothea Golze' via cp2k
Dear Ivan,
the q-TIP4P/F model of Habershon et al. contains an anharmonic
intramolecular O-H stretch, not a harmonic one (KIND G87). If I
remember correctly, in the original paper a quartic potential is described,
though there exist also an equivalent Morse-based parametrization.
Cheers,
Thomas
--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/aced09f1-216c-4430-ac1a-7aef56c06b3b%40googlegroups.com.
<test4.inp>
==============================
Thomas D. Kühne
Dynamics of Condensed Matter
Chair of Theoretical Chemistry
University of Paderborn
Warburger Str. 100
D-33098 Paderborn
Germany
Ivan Jackson
Mar 9, 2020, 5:17:50 PM3/9/20
to cp2k
Dear Prof. Kuehne,
I am really sorry for the confusion. I indeed changed the intramolecular O-H stretch to Quartic one as shown below:
&BOND
ATOMS OW HW
KIND QUARTIC
K 0.54192489 -0.9837788 0.9260195
R0 [angstrom] 0.9419
&END BOND
This is the corrected force field setup:
KIND QUARTIC
K 0.54192489 -0.9837788 0.9260195
R0 [angstrom] 0.9419
&END BOND
&BOND
ATOMS OW MW
KIND G87
K 0.0000000000000000E+00
R0 3.3070207325498757E-01
&END BOND
&BOND
ATOMS HW HW
KIND G87
K 0.0000000000000000E+00
R0 3.3070207325498757E-01
&END BOND
&BOND
ATOMS HW MW
KIND G87
K 0.0000000000000000E+00
R0 3.3070207325498757E-01
&END BOND
&BEND
ATOMS HW OW HW
KIND G87
K 0.139997992060672
THETA0 [deg] 107.40
I have attached the proper input file now.
Apologies for the confusion,
Best Regards,
Ivan Jackson
To unsubscribe from this group and stop receiving emails from it, send an email to cp...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/aced09f1-216c-4430-ac1a-7aef56c06b3b%40googlegroups.com.
<test4.inp>
Reply all
Reply to author
Forward
0 new messages