Which KSpace style is compatible with pair_style lj/charmm/coul/charmm?

LAMMPS LAMMPS Beginners
1

Just wonder which KSpace style is compatible with pair_style lj/charmm/coul/charmm ?
pair_style lj/charmm/coul/charmm 12.0 14.0 12.0 14.0

I tried

kspace_style ewald 1e-5
kspace_style pppm 1e-4
kspace_style msm 1e-4

All lead to errors like:
ERROR: KSpace style is incompatible with Pair style (src/kspace.cpp:199) ## if ewald
ERROR: KSpace style is incompatible with Pair style (src/kspace.cpp:201) ## if pppm
ERROR: KSpace style is incompatible with Pair style (src/kspace.cpp:203) ## if msm

2

Have you tried without any kspace_style command? Isn’t lj/charmm/coul/charmm dealing with electrostatic only in real space?

3

I see. Thanks!

4

None. This uses the (obsolete) smoothed cutoff coulomb.
Have you tried lj/charmm/coul/long (for old CHARMM force fields)? or lj/charmmfsw/coul/long (for current CHARMM force fields)?

Please make an effort to study the documentation more closely. It should be all explained there.

5

Hi Axel,
Thanks for your reply. I did read carefully the documentation at:
https://docs.lammps.org/pair_charmm.html#
On this webpage, there is no content about “real space”. And there are only two sentences about “kspace”:

  1. " All the styles with coul/long style are part of the KSPACE package."
  2. " Styles lj/charmm/coul/long and lj/charmm/coul/msm compute the same formulas as style lj/charmm/coul/charmm and style lj/charmmfsw/coul/long computes the same formulas as style lj/charmmfsw/coul/charmmfsh, except that an additional damping factor is applied to the Coulombic term, so it can be used in conjunction with the kspace_style command and its ewald or pppm or msm option."

The 2nd sentence is too long, and is really hard to understand. I suggest to split it into several short sentences, and mostly, I think some commas “,” are missing in the current sentence. Anyway, it doesn’t mention anything about the relationship about kspace and coul/charmm. How am I supposed to know it “uses the (obsolete) smoothed cutoff coulomb.” :slight_smile:

I searched “smoothed” on this webpage, and it only mentioned lj/charmm means that “the Lennard-Jones portion of the pair interaction is smoothed to 0.0 at the cutoff.” However, nothing is said about coul/charmm.

Actually, my real purpose is to reproduce the results of a coarse-grained (CG) force field (FF) in the format of the AMBER force field family and developed using the AMBER program. I searched the document of LAMMPS, and I only found the following:
https://docs.lammps.org/Howto_bioFF.html
“8.4.2. CHARMM, AMBER, COMPASS, DREIDING, and OPLS force fields.”
“CHARMM and AMBER”
Unfortunately, this page is not helpful for the purpose of applying AMBER force field in LAMMPS, except that it mentions “special_bonds amber” which is about the 1-4 scaling factor for lj and coul.

I found that the bonded energies (bonds, angles, dihedrals) can be exactly matched/re-produced by LAMMPS compared to AMBER by the following setting:
units real
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style charmm
(The old tool amber2lammps.py in the source code lammps-29Aug2024/tools/amber2lmp/ needs to be updated about the part of converting dihedral parameters)

However, the non-bonded energies (lj and coul) have significant discrepancies even I tried to use the same cut-off values. The density of pure organic liquids have a difference of about 4% under NPT ensemble at 298 K and 1 atm.

In order to debug, I build a simple system of 5 small CG molecules (5 particles/molecule, in total 25 CG atoms), in a cubic simulation box of 4nmx4nmx4nm, run 0 steps in both AMBER and LAMMPS, print out the coordinates (x, y, z) of each atom, the forces (fx, fy, fz) on each atom, and the bonded and non-bonded energy terms of the whole system. The printed coordinates are the same; the printed bonded energies for bonds, angles, and dihedrals are the same. However, the summed energies of VDW and COUL are not.
From AMBER, E(VDW) = VDWAALS + “1-4 NB”, E(COUL) = EELEC + “1-4 EEL”
From LAMMPS: E(VDW) = evdwl, E(COUL) = ecoul + elong

In AMBER, I set:
!! Nature and format of the output
ntxo = 1, !! ntox=1, formated ASCII for file “restrt”; ntxo=2, default, NetCDF file.
ioutfm = 0, !! ioutfm=0, formatted ASCII trajectory; ioutfm=1 (default) binary NetCDF trajectory
iwrap = 1, !! iwrap=1, coordinates to restart & trajectory file will be “wrapped”; default 0
ntpr = 1, !! frequency of printing energy information to files “mdout” and “mdinfo”, default 50
ntwx = 1, !! frequency of writing coordinate trajectory, default=0
ntwr = 1, !! frequency of writing “restrt” file, default=nstlim.
ntwf = 1 !! save the forces to a separate force trajectory file
!! Molecular dynamics
dt = 0.001, !! time step in ps
nstlim = 0, !! number of MD steps
!! Temperature regulation
ig = -1, !! the seed for the psedo-random number generator, affect starting velocity, Langevin dynamics or Andersen coupling.
!! ig=-1 (default), random seed based on current date and time, for ntt=2/3
ntt = 3, !! switch for temperature scaling
!! ntt=0, NVE ensemble
!! ntt=1, weak-coupling algorithm
!! ntt=2, Andersen temperature coupling scheme.
!! ntt=3, temperature, use Langevin dynamics with collision frequency gamma_ln
gamma_ln = 2.0, !! the collision frequency for thermostat (default 0). when ntt = 3, around 2 to 5 ps^-1
temp0 = 298.15, !! reference temperature (K), default 300
tempi = 298.15, !! initial temperature. default 0.0, velocities calculated from forces
!! When ntx<3 (for initial dynamics run), assign velocities from Maxwellian distribution.
!! When ntx>3, tempi has no effect.
tautp = 0.5, !! tautp, time constant (in ps) for heat bath coupling if ntt=1. default 1.0.
!! Pressure regulation
ntp = 1, !! flag for constant pressure dynamics
!! Note: TI is incompatible with ntp != 0, 1
!! ntp=0, default, no pressure scaling
!! ntp=1, pressure, isotropic position scaling
!! ntp=2, anisotropic (x-,y-,z-) pressure scaling
pres0 = 1.0, !! pres0=1.0 (default), reference pressure in units of bars
taup = 1.0, !! taup, default 1.0, pressure relaxation time in ps
!! Potential function parameters
ntb = 2, !! periodic boundary condition
!! ntb=0, no periodicity is applied and PME is off (default when igb>0)
!! ntb=1, constant volume (default when igb=0 & ntp=0)
!! ntb=2, constant pressure (default when ntp>1)
dielc = 1.0, !! dielectric multiplicative constant (default 1.0)
cut = 12.0, !! nonbonded cutoff in Angstroms (default 8.0 when igb==0)
igb = 0, !! flag for using implicit solvent models, default 0

In LAMMPS, I tried setting:
pair_style lj/charmm/coul/long 12.0 14.0 12.0
special_bonds amber
kspace_style ewald 1e-5
thermo 1
thermo_style custom step evdwl ecoul elong epair etail ebond eangle edihed eimp emol bonds angles dihedrals impropers
dump myForce all custom 1 myForce.dump id type x y z fx fy fz
run 0

If I change the pair_style to “pair_style lj/cut/coul/long 12.0 14.0” or “pair_style lj/charmmfsw/coul/long 12.0 14.0”, the discrepancies become even larger (run 0 step, compare the corresponding energies )

I found that If I want to get close values of E(VDW) and E(COUL) as “cut = 12.0” in AMBER, I need to set “pair_style lj/charmm/coul/long 14.0 20.0 13.12” in LAMMPS, so that the discrepancies of E(VDW) < 0.0001, E(COUL) < 0.0007, the discrepancies of forces (fx, fy, fz) on each atom < 0.005.

So, could you please kindly give some comments and suggestions on how to set up the pair_style (kspace), the ensemble, etc, in LAMMPS, in order to reproce the AMBER results as much as possible?

6

I cannot. I have not run anything with Amber in decades. So I honestly don’t know. Please also note that you will not be able to reproduce absolute energies as those are arbitrary and can have offsets, especially for the kspace parts, depending on the details of the implementation. You can only compare forces that are independent of them.

As you have seen, the support for force fields like CHARMM or Amber or GROMOS and alike is very shaky and often (very) outdated. Most of the tools and information goes back to the early days of LAMMPS, some 20 years ago and the the people with most experience in those have long ago left the LAMMPS development team, or were external LAMMPS contributors that just contributed tools some time ago and then moved on. Some contributors to LAMMPS even (sadly) passed on over the years.

I’ve been meaning to modernize all of this for a long time and set up test cases and explain how to convert inputs and parameters from one force field to the other and from one MD code to the other, but that is a lot of work and requires more time that I currently have available to work on LAMMPS. Particularly the CHARMM support is way outdated and does not take advantage of facilities that became available long after the implementation of the pair styles.
I do receive some financial support through a contract that Sandia has with Temple for my services, but that only pays for a fraction of my time, and that fraction is not even enough for managing current contributions and releases etc. through GitHub and responding to questions here to the best of my abilities (for as long as people asking questions are not trying to make me their drop-in adviser and ask me to tell them how to do their research).

Bottom line, under the current circumstances and with the fact that there is no revenue stream for LAMMPS and no sponsorship or similar available, the only way to get improvements, is to do them yourself and contribute them as a pull request. This is the price you and everybody else needs to pay these days for not having to pay any money to get access to all of LAMMPS.

Or put differently, if you absolutely want 100% Amber compatibility, you are better off using Amber’s pmemd instead of LAMMPS.

7

The good news is that you have achieved what you wanted, i.e. matched forces and energies between AMBER and LAMMPS, based on what you wrote above. Since these force fields involve taper functions at the VdW cutoff it’s not surprising that you have to use slightly different cutoffs to get exactly identical numerical behaviour between different implementations.

I would only suggest increasing the KSpace accuracy and seeing whether that has any effect on your results.

8

Axel,

As always, thanks for your long-time contributions to the communities of several computational programs and tools. ∠(‘-’ ) And thanks for your reply.

9

Hi srtee, thanks for your suggestion.

10

BTW: you should check out InterMol and ParmED for this purpose. Links are here: Pre/Post Processing Tools for use with LAMMPS

11

Axel, Thanks for the suggestion! I’ll check InterMol and ParmED.

12

all pair/angle/dihedral CHARMM styles and fix cmap were modernized with KOKKOS in January and November 2024. TLDR: generate input with https://charmm-gui.org/ and see examples/charmmfsw for PDB 1HVN. ive run CHARMM simulations with up to ~1M atoms on ~6000 cores for 7 days without issues.

ive been working on FitSNAP more lately but im still around. if a CHARMM question comes up please tag @alphataubio so i can help.

modern CHARMM KOKKOS pair/angle/dihedral styles have 66-89% unit test coverage which is higher than overall LAMMPS src at 38%.

13

Those changes do not address the fundamental issues with CHAMM support in LAMMPS that lies in the error-prone handling of 1-4 interactions and the fact that dihedrals should be using some variant of mult/harmonic instead of defining the same dihedral interaction multiple times with different parameters to recover the case of dihedrals with multiple multiplicities.

The stuff about the switching functions is an obsession with some folks in the CHARMM community and at the end something that in my personal opinion doesn’t matter given the limitations of a classical force field in general and the required underlying protocols for the force field parameterization. It was much more of a concern in the old days with (much) slower computers when people would reduce cutoffs to speed up calculations.

14

@akohlmey is dihedral_style fourier instead of dihedral_style multi/harmonic the correct way to resolve your concern ? or am i missing something subtle about the physics ??

ie. replace dihedral_style charmmfsw :

Dihedral Coeffs

           1       0.1000        1          180 1.00 # C CT1 CT1 CT2
           2       0.5200        2          180 0.00 # C CT1 CT1 CT2
           3       0.2900        3          180 0.00 # C CT1 CT1 CT2
[...]

with dihedral_style fourier :

dihedral_coeff 1  3  0.1 1 180  0.52 2 180  0.29 3 180 # C-CT1-CT1-CT2
15

The problem lies deeper in the implementation. With CHARMM compatible dihedral interactions, the 1-4 non-bonded terms may be fully excluded or scaled. The current implementation thus computes those non-bonded interactions as part of the dihedral style, not the pair style (where they belong).

Thus the CHAMM dihedral styles have an additional non-bonded weight factor as parameter. But this is not only to realize the different scaling of the 1-4 non-bonded interactions, but also to avoid computing the non-bonded contribution multiple times. This applies for example to rings with 6 or fewer members. For a 6-membered ring you have two dihedrals covering each pair of atoms at the opposite edge of the ring (1-4, 2-5, 3-6), so a weight of 0.5 needs to be applied to the full non-bonded interaction since you add it twice. For 5- or 4-membered rings, each 1-4 interaction is also a 1-3 or 1-2 interaction and those exclusions take precedence and thus the non-bonded part of the dihedral must be always fully excluded. Of course the weight factor of 0.0 is also needed if you specify the same dihedral a second or third time to describe dihedrals with more than one multiplicity, as they are present in the CHARMM parameter files. This can cause unexpected side effects through all kinds of functionality in LAMMPS.

This is mostly for historical reasons, since this functionality was added a long time ago for a now obsolete version of the CHARMM force field (19?), and when the neighbor list generation in LAMMPS was much simpler, but nobody has bothered to address this properly (instead they worry about that change in the switching function that happened since and the fact that CHARMM uses a slightly different conversion factor for Coulomb interactions (by about 0.0025%, while charges in the topology files are given with only two digits after the decimal point, so with an accuracy of about 1-2%)).

If I would have to do serious simulations with the CHARMM force field, I would stay away from LAMMPS and use either NAMD or Gromacs. They are both faster than LAMMPS, too, especially when running on CPUs.

16

what is the proper way to implement this ? what pair style in LAMMPS is the best example how to properly compute nonbonded interactions in the pair style where they belong ?

i see this code in dihedral_charmmfsw.cpp so how would this get moved to pair style properly:

    // 1-4 LJ and Coulomb interactions
    // tally energy/virial in pair, using newton_bond as newton flag

    if (weight[type] > 0.0) {
      itype = atomtype[i1];
      jtype = atomtype[i4];

      delx = x[i1][0] - x[i4][0];
      dely = x[i1][1] - x[i4][1];
      delz = x[i1][2] - x[i4][2];
      rsq = delx * delx + dely * dely + delz * delz;
      r2inv = 1.0 / rsq;
      r6inv = r2inv * r2inv * r2inv;

      // modifying coul and LJ force and energies to apply
      //   force_shift and force_switch as in CHARMM pairwise
      // LJ interactions between 1-4 atoms should usually be
      //   for r < cut_inner, so switching not applied

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

please explain these unexpected side effects and how to make sure to avoid avoid them in a new proper implementation in LAMMPS.

im willing to address this properly.

im not willing to do that. NAMD is not open source so thats a killer for me, and Gromacs has not updated unit testing in their repo since 8 years ago which is a huge red flag for me. i also looked at the gromacs code and it’s not even close to the software :germany: engineering quality in LAMMPS.

17

There is no precedent for this. This is part of the challenge. Plus having to overcome the “why change something that is already working?” barrier.

I am sorry, but I have to stop here. I see this is heading down a path that I want to avoid for personal reasons.