pSPICA
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 .
Note that
to use the LJ12-5 function form, LAMMPS patch 3 AUG 2022 or later version
is required.
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
wat2polar
command of
spica-tools.
The command line will be:
$ cg_spica wat2polar spica_conf.pdb pspica_conf.pdb
where the first and second argument following wat2polar 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_lmp command for the initial configuration.
The script requires parameter and topology files of the pSPICA force field.
For this system, we prepared
PDMPC.top and
PWAT.top
as the topology files,
and
pSPICA_lipid.prm 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 = 80).
setup_lmp
will provide 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
file.
$ cg_spica setup_lmp PDMPC.top 128 PWAT.top 1152 pSPICA_lipid.prm pspica_conf.pdb
Check the generated two LAMMPS input files,
DATA.FILE and
PARM.FILE.
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,
fix shake.
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
...
The
b 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 PARM.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,
PARM.FILE, 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 the "wat2polar" command to replace the uncharged water of the SPICA FF in the initial configuration,
system.cg.pdb,
with the polar water of pSPICA:
$ cg_spica wat2polar system.cg.pdb p_system.cg.pdb
Then we will excute "setup_lmp" 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.
The
ENM command used in the SPICA protein tutorial
has an option,
-pspica, to assign partial charges for pSPICA. So with
protein.cg.pdb, you can generate the
protein topology file where values of partial charges are 0.5990 using the following command:
$ cg_spica ENM protein.cg.pdb p_protein.cg.top -pspica
p_protein.cg.top is a protein topology file for the pSPICA FF.
Then we will use "setup_lmp" 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).
$ cg_spica setup_lmp p_protein.cg.top 1 PDOPC.top 128 PWAT.top 2134 PCLA.top 4 pSPICA_protein.prm p_system.cg.pdb
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, PARM.FILE, and the LAMMPS setup file in.npt.
Make sure the command "fix shake" are properly set in the input, mentioned
above.