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”:
- " All the styles with coul/long style are part of the KSPACE package."
- " 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.” 
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?