SPC/E in force_fields/spce_oplsaa.lt

[feifzhou]

The spc_oplsaa.lt file and the one in the example are inconsistent.

diff moltemplate/moltemplate/force_fields/spc_oplsaa.lt moltemplate/examples/all_atom/force_field_OPLSAA/waterSPC_using_OPLSAA/moltemplate_files/spc.lt

Notably, the constraint is missing in the latter.

1c1
< # file "spc_oplsaa.lt" 
---
> # file "spce.lt" 
6,8d5
< #
< # This is the version of the SPC water model suitable for use in a simulation
< # using the OPLSAA force field (as implemented in the "oplsaa.lt" file).
10d6
< import "oplsaa.lt"  # <-- defines OPLSAA, @atom:76 and @atom:77
12c8,10
< SPC_oplsaa inherits OPLSAA {
---
> import "oplsaa.lt"
> 
> SPC inherits OPLSAA {
14,16c12,14
<   # Atom types from "oplsa.lt"
<   # @atom:76 <--> "SPC Water O"
<   # @atom:77 <--> "SPC Water H"
---
>   # Atom types from "oplsaa_lt_generator/oplsaa_subset.prm"
>   # @atom:76  <-->  OW  "SPC Water O"
>   # @atom:77  <-->  HW  "SPC Water H"
19,21c17,19
<     $atom:O  $mol:. @atom:76 -0.82   0.0000000 0.00000 0.000000
<     $atom:H1 $mol:. @atom:77  0.41   0.8164904 0.00000 0.5773590
<     $atom:H2 $mol:. @atom:77  0.41  -0.8164904 0.00000 0.5773590
---
>     $atom:O  $mol:. @atom:76 -0.8200 0.0000000 0.0000 0.000000
>     $atom:H1 $mol:. @atom:77 0.4100 0.8164904 0.0000 0.577359
>     $atom:H2 $mol:. @atom:77 0.4100 -0.8164904 0.0000 0.577359
25,26c23,24
<     $bond:OH1  $atom:O $atom:H1
<     $bond:OH2  $atom:O $atom:H2
---
>     $bond:OH1 $atom:O $atom:H1
>     $bond:OH2 $atom:O $atom:H2
29,39c27,35
<   write_once("In Constraints") {
<     # Define a group for the spc water molecules:
<     group spc type  @atom:76  @atom:77
< 
<     # Constrain the angles and distances in spc water:
<     fix fRattleSPCE spc rattle 0.0001 10 100 b @bond:OH a @angle:HOH
< 
<     # Remember to put this command in your LAMMPS input script:
<     #  include system.in.constraints
<     # ...AFTER minimization and after all integration fixes.
<   }
---
> } # end of definition of "SPC" water molecule type
> 
> 
> 
> 
> 
> 
> 
> 
41d36
< } # end of definition of "SPCE" water molecule type
42a38,95
> 
> ###################### old version (SPCE) ######################
> #
> #SPCE {
> #
> #  write("Data Atoms") {
> #    $atom:O  $mol:. @atom:O -0.8476 0.0000000 0.0000 0.000000
> #    $atom:H1 $mol:. @atom:H 0.4238 0.8164904 0.0000 0.577359
> #    $atom:H2 $mol:. @atom:H 0.4238 -0.8164904 0.0000 0.577359
> #  }
> #
> #  write_once("Data Masses") {
> #    @atom:O 15.9994
> #    @atom:H 1.008
> #  }
> #
> #  write("Data Bonds") {
> #    $bond:OH1 @bond:OH $atom:O $atom:H1
> #    $bond:OH2 @bond:OH $atom:O $atom:H2
> #  }
> #
> #  write("Data Angles") {
> #    $angle:HOH @angle:HOH $atom:H1 $atom:O $atom:H2
> #  }
> #
> #  write_once("In Settings") {
> #    bond_coeff   @bond:OH     harmonic               1000.0  1.0 
> #    angle_coeff  @angle:HOH   harmonic               1000.0  109.47
> #    pair_coeff   @atom:O @atom:O  lj/cut/coul/long   0.1553  3.166 
> #    pair_coeff   @atom:H @atom:H  lj/cut/coul/long   0.0     2.058
> #    # Note: Since we are using RATTLE constraints, the bond and angle strength
> #    # parameters ("600.0", "75.0") do not matter. But the equilibrium bond
> #    # length ("1.0") and equilibrium angle ("109.47") does matter. LAMMPS
> #    # obtains these numbers from the bond_coeff and angle_coeff commands above.
> #  }
> #
> #  write_once("In Constraints") {
> #    group spce type  @atom:O  @atom:H
> #    fix fRattleSPCE spce rattle 0.0001 10 100 b @bond:OH a @angle:HOH
> #    # Remember to put this command in your LAMMPS input script:
> #    #  include system.in.constraints
> #    # ...AFTER minimization and after all integration fixes.
> #  }
> #
> #  write_once("In Init") {
> #    # -- Default styles (for solo "SPCE" water) --
> #    units        real
> #    atom_style   full
> #    # (Hybrid force fields were not necessary but are used for portability.)
> #    pair_style   hybrid lj/cut/coul/long 10.0
> #    bond_style   hybrid harmonic
> #    angle_style  hybrid harmonic
> #    kspace_style pppm 0.0001
> #    pair_modify  shift yes
> #  }
> #
> #} SPCE
> ###################################################################

[feifzhou]

P.S. I had a problem using moltemplate/moltemplate/force_fields/spc_oplsaa.lt

$ cat system.lt 
import "spc_oplsaa.lt"     # <- defines the "SPC" (water) molecule type (uses OPLSAA)

write_once("Data Boundary") {
   0.0  41.40  xlo xhi
   0.0  41.40  ylo yhi
   0.0  41.40  zlo zhi
}

waters = new SPC_oplsaa [12].move(0.00, 0.00, 3.45) 
                 [12].move(0.00, 3.45, 0.00) 
                 [12].move(3.45, 0.00, 0.00)

Now there is an error

$ moltemplate.sh system.lt 
moltemplate.sh v2.20.14 2022-9-09

lttree_check.py v0.81.2 2021-5-24
########################################################
##            WARNING: atom_style unspecified         ##
## --> "Data Atoms" column data has an unknown format ##
##              Assuming atom_style = "full"          ##
########################################################
lttree_check.py:    parsing the class definitions... done
lttree_check.py:    looking up classes... done
lttree_check.py:    looking up @variables... done

---------------------------------------------------------------------
     Syntax error: Missing bond coeff.

  No coeffs for the "@/bond:SPC_oplsaa/OH" bond type have been
defined, but a reference to that bond type was discovered
near "spc_oplsaa.lt", line 34.   Check this file and also check
your "bond_coeff" commands or your "Data Bond Coeffs" section.
---------------------------------------------------------------------
(To continue anyway, run moltemplate using the "-nocheck" argument.)

[jewettaij]

Thank you Fei Zhou for catching this mistake. I have updated the "spc_oplsaa.lt" and "tip3p_oplsaa.lt" files. I posted the commit today. It's embarrassing. I never tested either of these files and they both had bugs. Hopefully they work better now. I also added some comments to both files to better explain how they should be used.

Also: The SPC and SPC/E (also known as "SPCE") water models are not the same. They have different parameters. If you are curious, I think both water models are explained in this paper: Berendsen, Grigera, Straatsma, J Phys Chem, 91, 6269-6271 (1987)

[jewettaij]

These changes were made in the commit I posted today (e2e8782).

[jewettaij]

I will close this issue now. Feel free to reopen it.