Here comes a few key pointer how to setup a hybrid AA/CG simulation with Lammps. For specific examples, please consult the various tutorials.
The atom style
The atom style should be a hybrid of three styles to get the molecular id, charge, dipole as well as density and diameter properties.
atom_style hybrid sphere dipole molecular
The pair style and pair coefficients
Here it starts to be a little bit complicated, because we want to use an rRESPA integrator and such we need to have a pair style for CG-AA, CG-CG and AA-AA forces. In addition, you might have different AA force fields and a such need to separate these as well.
The pair style of the CG-CG and CG-AA interactions are lj_sf_dipole_sf. For AA-AA, it is mostly tested with an Amber/Charmm type of force field, which uses the style lj/charmm/coul/long. (For some molecules you might want to use the lj/charmm/coul/long/14, to allow for individual 1-4 scaling).
Putting everything together, this becomes
pair_style hybrid lj/sf/dipole/sf 12.0 lj/sf/dipole/sf 12.0 lj/charmm/coul/long 11.0 12.0
where 12.0 and 12.0 11.0 are the non-bonded cut-offs in Å.
Then we need to consider the scaling of the 1-2, 1-3 and 1-4 interactions. ELBA is parameterised with a scaling of 1.0 for 1-3 and 1-4 interactions, whereas Amber uses 0.0 for 1-3 and 0.5/0.83 for 1-4. This can be accomplished with these commands
special_bonds lj/coul 0.0 0.999999 0.999999
pair_modify pair lj/charmm/coul/long special lj 0.0 0.0 0.5
pair_modify pair lj/charmm/coul/long special coul 0.0 0.0 0.8333
The first command sets the ELBA parameters and the other commands sets the Amber scaling for only the lj/charmm/coul/long/ style.
The Amber pair style requires a long-range solver, i.e. we will use a particle-particle particle-mesh Ewald solver
kspace_style pppm/cg 1.0e-5
We then need to specify the pair parameter, i.e. Lennard-Jones parameters (and additional electrostatic scaling) with the pair_coeff commands. These commands, which will not change from simulations to simulations, are best put in a separate “include”-file. They are usually generated by scripts, because it will be very tedious to add all the parameters by hand.
An example is
pair_coeff 1 8 lj/sf/dipole/sf 1 0.30578 3.15000
where 1 8 are the atom types. lj/sf/dipole/sf is the pair style for this pair and the 1 immediately after indicates that this should interaction should be computed with the first style in the hybrid command. The last two numbers are the Lennard-Jones epsilon and sigma parameters.
For CG-CG interactions an example is
pair_coeff 7 7 lj/sf/dipole/sf 2 0.83700 4.50000
where the 2 immediately following the pair style indicates that this interaction should be computed by the second style in the hybrid command. This will become important when we look at the rRESPA settings below.
For some lipid beads, the charge-dipole interaction needs to be scaled in order to recover proper behavior, this is accomplished with a line like this
pair_coeff 5 20 lj/sf/dipole/sf 1 0.32340 3.99984 scale 1.750
Other styles
The styles for the bonds, angles and dihedrals need to look something like this
bond_style harmonic
angle_style hybrid cosine/squared dipole harmonic
dihedral_style harmonic
The ELBA angle style is a hybrid of cosine/squared and dipole.
If you mix two AA force fields, you might need to use a hybrid style for the dihedrals.
dihedral_style hybrid charmm harmonic
Setting up the integrator
For the hybrid simulations it is recommended to use the rRESPA integrator, which allows the CG-CG pair forces to be propagated at much longer time steps than the other forces.
First you set the time step for the CG-CG pair forces
timestep 8.0
and then you set the run style, e.g.
run_style respa 2 4 hybrid 1 2 1 kspace 1
This tells LAMMPS to use rRESPA with 2 levels. The time step of the outer layer is set by the timestep command and the time step of the inner level is set as a factor of the outer step. Here, we use a factor of 4, which means that the inner forces are computed 4 times for every computation of the outer forces, i.e. the inner time step is 2 fs.
The hybrid keyword is followed by the level at which each of the pair styles specified above will be computed. Here we set that the first instance of lj/sf/dipole/sf and lj/charmm/coul/long is propagated with a 2 fs time step, and that the second instance of lj/sf/dipole/sf is propagated with a 8 fs time step. As we specified with the pair_coeff commands that CG-CG interactions are computed with the second instance, we have the desired effect of propagating them with a longer time step. (Naturally, we could have reversed the numbers in the pair_coeff and run_style commands.)
Setting up fixes
First, we will set up different groups that differentiate the AA and CG particles. This can be done by atom type for instance
group cg type 1:7
group aa type > 7
as there are 7 ELBA bead types.
Then to propagate the positions and velocities of the particles, we need to setup two fixes, one for the AA atoms and one for the CG beads.
fix integrateSolute aa nve
fix integrateCG cg nve/sphere update dipole/dlm
Here, we have used the nve integrators, which means we also need to setup thermostats and barostat. Alternatively, we can use the npt and npt/sphere fixes.
Load balancing
To improve the parallell efficiency of the simulations, it is recommended to use an improved load balancing, weighting the AA particles more (as they are propagated more often).
balance 1.1 shift xyz 5 1.1 weight group 1 aa 4.0
do this at the beginning of the run, but then does not update it.
fix loadbalance all balance 100 1.1 shift xyz 5 1.1 weight group 1 aa 4.0
on the other hand updates the balance every 100 step. You should benchmark this and the load factor of 4, before you run any large simulations.