pSPICA is a CG FF designed with a polar CG water model.
The CG modeling method of pSPICA is the same as the SPICA FF.
The LJ12-5 function is adopted to reproduce the experimental density, surface tension,
and dielectric permittivity of water.
A brief summary of the model is given here
and the publication with full details here
For simulations with LAMMPS there is already a script to use the LJ12-5 function
by changing the LAMMPS parameter file generated by setup_lammps
The tutorial for generating LAMMPS parameter files can be found here
There is another script to change the uncharged CG water of the SPICA FF to the polar CG water of pSPICA.
Example for lipid membrane system:
Here we will show you how to use the pSPICA to simulate a lipid membrane system
composed of 128 DMPC lipids and 1152 polar CG water particles with LAMMPS,
and start the simulation from an initial configuration obtained through the SPICA tutorial
for a lipid membrane system.
The equilibrated configuration, spica_conf.pdb
, has the standard SPICA CG water, representing three water molecules as
a single CG particle.
So, we replace the SPICA CG water in the configuration with the pSPICA polar CG water using
The command to excute the code is:
$ python WAT2PWAT.py spica_conf.pdb pspica_conf.pdb
where the first argument is the python code. The second and third arguments are input (using the SPICA water) and output (using the pSPICA water, pspica_conf.pdb
) configuration files, respectively.
Next, to generate LAMMPS input files, we run the setup_lammps
script for the initial configuration. The script requires parameter and
topology files of the pSPICA force field. For this system, we prepared
as the topology files,
as the parameter file of pSPICA lipids and water.
We put the prefix of pSPICA topology files, "P"
, to distinguish those from SPICA ones. Please be careful the magnitude
of partial charge in pSPICA is different from that in SPICA. For example,
the CG segment PH, corresponding to phosphate group, has a partial charge
of -0.5990 in pSPICA topology files. On the other hand, the segment has
a partial charge of +0.1118 in SPICA topology files. These values are taken
to consider the relative permittivity implicitly (pSPICA:εr
= 3.2 and SPICA:εr
"setup_lammps" can be easily installed by the following commands:
$ tar xvf setup_lammps_v2.tar
$ cd setup_lammps_v2
Then, setup_lammps will produce the LAMMPS readable files from topology,
parameter, and configuration files. The number of molecules contained in
the target system must follow topology file names in the command. In addition,
the topology file names must be in sequence indentical to the system configuration
$ setup_lammps PDMPC.top 128 PWAT.top 1152 pSPICA_lipid.prm pspica_conf.pdb
Check the generated two LAMMPS input files, DATA.FILE
and PARM.FILE. After that, we modify the obtained PARM.FILE to use tabulated
potential with the pair_style table
The LJ12-5 function, one of the interactions between CG water molecules, cannot be adopted
currently through the pair_style lj/sdk
$ convert_PARM.py PARM.FILE pPARM.FILE
The generated pPARM.FILE
should look like:
# Generated by setup_lammps
pair_style hybrid/overlay lj/sdk 15.0 coul/long 15.0 table linear 2999
angle_style hybrid sdk harmonic
pair_coeff 1 1 lj/sdk lj9_6 0.5650 5.4500 # NC NC
pair_coeff 1 1 coul/long # NC NC
pair_coeff 8 8 table 12_5ww.t TABLE 15.0 # WO WO
pair_coeff 8 8 coul/long # WO WO
The pair_style hybrid/overlay
enables us to use multiple pair styles. Using the style, we accomplish to describe LJ interactions, LJ12-4, LJ9-6 and LJ12-5,
with the combination of lj/sdk and table. 12_5ww.t
is a file containing tabulated energy and force values of the LJ12-5 interaction between CG water. Interaction among CG water and ion in pSPICA is described by LJ12-5. Although in this tutorial we are not going to use CG ions (CG sodium and chloride), their topology and tabulated potential files are prepared here
for your information.
The .tar file includes the following files :
- XXX.top : topology file of CG water (XXX = PWAT : CG water, PSOD : CG sodium, PCLA : CG chloride)
- 12_5xy.t : table file to apply tabulated potential (x, y = w : water, s : sodium, or c : chloride)
The remaining file to perform the simulation with pSPICA is a LAMMPS setup file.
We use in.npt
as the setup file
to optimize the configuration and simulate the system with the NPT ensemble.
Please make sure that the file includes a mandatory command,
pSPICA CG water represented by two particles having oposite charge each other has
a rigid bond to connect the particles.
The two commands in the in.npt should look like:
fix Fix0 all shake 1.0e-6 25 0 b 9
of the "fix shake" command lists the bond types that will be constrained.
We have to be careful about the last argument of the command
because it corresponds to the bond index that is fixed.
The bond index and parameter in the pPARM.FILE look like:
bond_coeff 7 8.4000 3.4800 # GL EST2
bond_coeff 8 8.4000 3.4800 # EST2 CM
bond_coeff 9 8.4000 3.4800 # WH WO
angle_coeff 1 sdk 3.1000 112.0000 lj9_6 0.5550 4.6200 # NC PH GL
We can see that the index of bond between WH-WO, composing the CG water, is 9 there.
This is why the last argument of "fix shake" is set to 9 in the in.npt.
Whenever you build new system configuration files and prepare input files
to carry out simulations with pSPICA, you have to check the bond index
of WO-WH (and CG ions), and change the last argument of "fix shake"
to the bond index number
. In addition, please be careful to write the "fix shake" command
between the "minimize" and "fix npt" command. Otherwise,
some errors in shake determinant or energy calculation will occur in the
begining of the simulations.
The LAMMPS simulation can now be run in the directory where the DATA.FILE,
pPARM.FILE, 12_5ww.t, and the in.npt are, with the following command:
$ lmp_mpi -in in.npt -log output_npt.log
Example for lipid-protein system:
We can also simulate lipid-protein systems using the pSPICA FF.
After gaining the initial configuration of a target lipid-protein system by following
tutorial for protein systems
we will apply "WAT2PWAT.py" to replace the uncharged water of the SPICA FF in the initial configuration,
with the polar water of pSPICA:
$ WAT2PWAT.py system.cg.pdb p_system.cg.pdb
Then we will excute "setup_lammps" using topology files for the pSPICA FF.
Although topology files for protein can be generated by following tutorial for protein systems
one must be careful about values of the charge in the pSPICA FF.
We need to change a partial charge value of 0.1118 to 0.5990 in CG protein topology files generated with "gen_top_elastic_network.py".
In addition, if the terminal backbone particles, "GBT", in the protein topology file have positive or negative charge,
we need to change the CG type "GBT" to the CG type "GBTP" (in case of positive charge) or "GBTN" (in case of negative charge).
For charge values in protein topology files, they can be easily modified like this:
$ sed -e "s/0.1118/0.5990/g" protein.cg.top > p_protein.cg.top
is a protein topology file for the SPICA FF.
is a protein topology file for the pSPICA FF.
Then we will use "setup_lammps" command to generate LAMMPS input files, similarly to tutorial for protein systems
Download the pSPICA protein parameter set (pSPICA_protein.prm)
and topology files of DOPC(PDOPC.top)
and CLA (PCLA.top)
$ setup_lammps p_protein.cg.top 1 PDOPC.top 128 PWAT.top 2134 PCLA.top 4 pSPICA_protein.prm p_system.cg.pdb
After that, we modify the obtained PARM.FILE, same as tutorial for
lipid membrane systems
$ convert_PARM.py PARM.FILE pPARM.FILE
are required to apply tabulated potentials between CG water and chloride, and among CG chloride, repectively.
A LAMMPS input in.npt
for the lipid-protein system is almost same as that for the lipid membrane system.
Then we can perform LAMMPS simulations of the lipid-protein system using
DATA.FILE, pPARM.FILE, 12_5ww.t, 12_5wc.t, 12_5cc.t, and the LAMMPS setup file in.npt.
Make sure the command "fix shake" are properly set in the input, mentioned above