shirtsgroup / InterMol

Conversion tool for molecular simulations
MIT License
186 stars 53 forks source link

Pairtypes are missing in the generated lammps input #314

Open maaren opened 7 years ago

maaren commented 7 years ago

When I run intermol with https://github.com/dspoel/liquid-simulations/blob/master/FF/GROMACS/CGenFF/TOP/GAS/1-butanol.top as input I get:

            type     input(gromacs)   output (lammps)     diff (lammps)
           Angle                n/a        1.42339730               n/a
   All dihedrals         4.20188051        4.20188032       -0.00000019
      All angles         1.42339729        1.42339730        0.00000001
      Coulomb-14       -10.78972151               n/a               n/a
       Potential        24.55184221       24.80523228        0.25339007
   Electrostatic        20.45497800       20.45497811        0.00000010
           LJ-14        -0.10460745               n/a               n/a
            Bond         0.38551519        0.38551519        0.00000000

when I remove the [ pairtypes ] section in 1-butanol.top and rerun I get:

            type     input(gromacs)   output (lammps)     diff (lammps)
           Angle                n/a        1.42339730               n/a
   All dihedrals         4.20188051        4.20188032       -0.00000019
      All angles         1.42339729        1.42339730        0.00000001
      Coulomb-14       -10.78972151               n/a               n/a
       Potential        27.70885244       27.70885254        0.00000011
   Electrostatic        20.45497800       20.45497811        0.00000010
           LJ-14         3.05240278               n/a               n/a
            Bond         0.38551519        0.38551519        0.00000000

So something goes wrong with 1-4 input of lammps. After reading: http://lammps.sandia.gov/doc/pair_charmm.html

I tried to use: pair_style lj/charmm/coul/charmm 19.999 20.0

and experimented with the pair_coeff parameters to get epsilon_14 and sigma_14 in there, but until now I have not got something that reproduces the energies with lammps.

maaren commented 7 years ago

I don't think it is very clearly documented, but I think that the LJ14 interactions are added in the dihedral potential. See: lammps-17Nov16/src/dihedral_charmm.cpp

  if (implicit) forcecoul = qqrd2e * q[i1]*q[i4]*r2inv;
  else forcecoul = qqrd2e * q[i1]*q[i4]*sqrt(r2inv);
  forcelj = r6inv * (lj14_1[itype][jtype]*r6inv - lj14_2[itype][jtype]);
  fpair = weight[type] * (forcelj+forcecoul)*r2inv;

  if (eflag) {
    ecoul = weight[type] * forcecoul;
    evdwl = r6inv * (lj14_3[itype][jtype]*r6inv - lj14_4[itype][jtype]);
    evdwl *= weight[type];
  }

I have managed to create a lammps charmm dihedral topology which reproduces the gmx dihedral energy for 1-butanol, but I still have not been able to get the lammps E_vdwl equal to gmx LJ 14+SR