git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11956 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2014-05-09 15:35:51 +00:00
parent db78f042ef
commit 944b898d71
64 changed files with 108920 additions and 0 deletions

View File

@ -0,0 +1,87 @@
Ethane to Methanol in Water
===========================
Example calculation of the difference in free energy of hydration upon
transforming ethane into methanol with LAMMPS using *compute fep* and
*fix adapt/fep*.
Ethane and methanol are represented by the OPLS-AA force field (1
molecule). Water is represented by the 3-site SPC/E model (360
molecules).
The strategy used to perform the alchemical transformation is the
following:
* The dual topology approach is used, therefore all the atoms of
ethane and methanol are present throughout the simulation, only some
of them are dummy sites at the endpoints of the
transformation. Masses and intramolecular terms (bond lengths,
angles, dihedrals) are not changed.
* Interactions of sites that are being created (from dummy sites) or
deleted (to become dummy sites) are treated using soft-core verions
of the Lennard-Jones and Coulomb potentials (*pair
lj/cut/coul/long/soft*) in order to avoid singularities. The
exponent of the coupling parameter lambda in the soft-core pair
potentials was in this example n = 1.
* Eletrostatic charges that are modified are varied linearly from the
initial to the final values. This keeps the overall charge of the
molecule constant, which is good for the long range electrostatics
(the coupling parameter lambda has no effect on the kspace terms).
The following directories contain input files and results for
calculations using free-energy perturbation (FEP), thermodynamic
integration (TI/FDTI) and Bennet's acceptance ratio methods:
* `mols` -- Molcule description files and force field database used to
create the initial configurations for the simulations `data.0.lmp`
and `data.1.lmp`
* `fep01` -- Calculation using FEP, multi-stage transformation of an
ethane molecule into methanol. Results in `fep01.lmp`
* `fep10` -- Calculation using FEP, multi-stage transformation of a
methanol molecule into ethane. Results in `fep10.lmp`
* `fdti01` -- Calculation using FDTI, transformation of an
ethane molecule into methanol. Results in `fdti01.lmp`
* `fdti10` -- Calculation using FDTI, transformation of a
methanol molecule into ethane. Results in `fdti10.lmp`
* `bar01` -- Calculation using BAR, 1-step transformation of an
ethane molecule into methanol. Results in `bar01.lmp`
* `bar10` -- Calculation using BAR, 1-step transformation of a
methanol molecule into ethane. Results in `bar10.lmp`
The free-energy profiles can be observed by plotting the values in the
third column of the results files. The Python scripts `fep.py`,
`nti.py`, `fdti.py`, and `bar.py` found in the `tools` directory can
be used to calculate the free-energy differences corresponding to the
above transformations:
fep.py 300 < fep01.lmp
fep.py 300 < fep10.lmp
nti.py 300 0.002 < fdti01.lmp
nti.py 300 0.002 < fdti10.lmp
fdti.py 300 0.002 < fdti01.lmp
fdti.py 300 0.002 < fdti10.lmp
bar.py 300 bar01.lmp bar10.lmp
The outputs are in kcal/mol and can be compared with the experimental
value of -6.93 kcal/mol and with simulation value from the literature:
-6.7 kcal/mol [Jorgensen, Ravimohan, J Chem Phys 83 (1985) 3050], -6.8
kcal/mol [Goette, Grubmüller, J Comp Chem 30 (2007) 447].
These example calculations are for tutorial purposes only. The results
may not be of research quality (not enough sampling, size of the step
in lambda or of the delta for numerical derivative not optimized, no
evaluation of ideal-gas contributions, etc.)

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,120 @@
# created by fftool
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
special_bonds lj/coul 0.0 0.0 0.5
read_data data.0.lmp
# read_restart restart.*.lmp
# reset_timestep 0
pair_style hybrid lj/cut/coul/long 10.0 10.0 &
lj/cut/coul/long/soft 1 0.5 10.0 10.0 10.0
pair_modify tail no
kspace_style pppm 1.0e-4
pair_coeff 1 1 none # C-CO C-CO
pair_coeff 1 2 none # C-CO C-D
pair_coeff 1 3 none # C-CO H-D
pair_coeff 1 4 none # C-CO H-H1
pair_coeff 1 5 none # C-CO OHD-
pair_coeff 1 6 none # C-CO HOD-
pair_coeff 1 7 lj/cut/coul/long 0.0000 1.0000 # C-CO Hw
pair_coeff 1 8 lj/cut/coul/long 0.1013 3.3286 # C-CO Ow
pair_coeff 2 2 none # C-D C-D
pair_coeff 2 3 none # C-D H-D
pair_coeff 2 4 none # C-D H-H1
pair_coeff 2 5 none # C-D OHD-
pair_coeff 2 6 none # C-D HOD-
pair_coeff 2 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # C-D Hw
pair_coeff 2 8 lj/cut/coul/long/soft 0.1013 3.3286 1.0 # C-D Ow
pair_coeff 3 3 none # H-D H-D
pair_coeff 3 4 lj/cut/coul/long/soft 0.0300 2.5000 1.0 # H-D H-H1
pair_coeff 3 5 none # H-D OHD-
pair_coeff 3 6 none # H-D HOD-
pair_coeff 3 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # H-D Hw
pair_coeff 3 8 lj/cut/coul/long/soft 0.0683 2.8131 1.0 # H-D Ow
pair_coeff 4 4 none # H-H1 H-H1
pair_coeff 4 5 none # H-H1 OHD-
pair_coeff 4 6 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # H-H1 HOD-
pair_coeff 4 7 lj/cut/coul/long 0.0000 1.0000 # H-H1 Hw
pair_coeff 4 8 lj/cut/coul/long 0.0683 2.8131 # H-H1 Ow
pair_coeff 5 5 none # OHD- OHD-
pair_coeff 5 6 none # OHD- HOD-
pair_coeff 5 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # OHD- Hw
pair_coeff 5 8 lj/cut/coul/long/soft 0.1625 3.1427 0.0 # OHD- Ow
pair_coeff 6 6 none # HOD- HOD-
pair_coeff 6 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # HOD- Hw
pair_coeff 6 8 lj/cut/coul/long/soft 0.0000 1.7792 0.0 # HOD- Ow
pair_coeff 7 7 lj/cut/coul/long 0.0000 1.0000 # Hw Hw
pair_coeff 7 8 lj/cut/coul/long 0.0000 1.0000 # Hw Ow
pair_coeff 8 8 lj/cut/coul/long 0.1554 3.1655 # Ow Ow
variable nsteps equal 500000
variable nprint equal ${nsteps}/500
variable ndump equal ${nsteps}/100
# variable nrestart equal ${nsteps}/10
variable temp equal 300.0
variable press equal 1.0
neighbor 2.0 bin
timestep 2.0
# velocity all create ${temp} 12345
thermo_style multi
thermo ${nprint}
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 4 5 a 6
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 500
run 100000
reset_timestep 0
variable dlambda equal 1.0
variable minusdl equal -1.0
variable dq1 equal 0.145*v_dlambda+0.180*v_dlambda
variable dq2 equal 0.180*v_dlambda
variable dq3 equal -0.060*v_dlambda
variable dq4 equal 0.040*v_dlambda-0.060*v_dlambda
variable dq5 equal -0.683*v_dlambda
variable dq6 equal 0.418*v_dlambda
compute cFEP all fep ${temp} &
pair lj/cut/coul/long/soft lambda 2*3 7*8 v_minusdl &
pair lj/cut/coul/long/soft lambda 3 4 v_minusdl &
pair lj/cut/coul/long/soft lambda 4 6 v_dlambda &
pair lj/cut/coul/long/soft lambda 5*6 7*8 v_dlambda &
atom charge 1 v_dq1 &
atom charge 2 v_dq2 &
atom charge 3 v_dq3 &
atom charge 4 v_dq4 &
atom charge 5 v_dq5 &
atom charge 6 v_dq6
fix fFEP all ave/time 1 1 100 c_cFEP[1] c_cFEP[2] file bar01.lmp
# compute cRDF all rdf 100 1 1
# fix fRDF all ave/time 20 100 ${nsteps} c_cRDF file rdf.lammps mode vector
# compute cMSD all msd
# fix fMSD all ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C C H H O H H O
# restart ${nrestart} restart.*.lmp
run ${nsteps}
write_restart restart.*.lmp
write_data data.*.lmp

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,120 @@
# created by fftool
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
special_bonds lj/coul 0.0 0.0 0.5
read_data data.1.lmp
# read_restart restart.*.lmp
# reset_timestep 0
pair_style hybrid lj/cut/coul/long 10.0 10.0 &
lj/cut/coul/long/soft 1 0.5 10.0 10.0 10.0
pair_modify tail no
kspace_style pppm 1.0e-4
pair_coeff 1 1 none # CO-C CO-C
pair_coeff 1 2 none # CO-C CD-
pair_coeff 1 3 none # CO-C HD-
pair_coeff 1 4 none # CO-C H1-H
pair_coeff 1 5 none # CO-C OH-D
pair_coeff 1 6 none # CO-C HO-D
pair_coeff 1 7 lj/cut/coul/long 0.0000 1.0000 # CO-C Hw
pair_coeff 1 8 lj/cut/coul/long 0.1013 3.3286 # CO-C Ow
pair_coeff 2 2 none # CD- CD-
pair_coeff 2 3 none # CD- HD-
pair_coeff 2 4 none # CD- H1-H
pair_coeff 2 5 none # CD- OH-D
pair_coeff 2 6 none # CD- HO-D
pair_coeff 2 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # CD- Hw
pair_coeff 2 8 lj/cut/coul/long/soft 0.1013 3.3286 0.0 # CD- Ow
pair_coeff 3 3 none # HD- HD-
pair_coeff 3 4 lj/cut/coul/long/soft 0.0300 2.5000 0.0 # HD- H1-H
pair_coeff 3 5 none # HD- OH-D
pair_coeff 3 6 none # HD- HO-D
pair_coeff 3 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # HD- Hw
pair_coeff 3 8 lj/cut/coul/long/soft 0.0683 2.8131 0.0 # HD- Ow
pair_coeff 4 4 none # H1-H H1-H
pair_coeff 4 5 none # H1-H OH-D
pair_coeff 4 6 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # H1-H HO-D
pair_coeff 4 7 lj/cut/coul/long 0.0000 1.0000 # H1-H Hw
pair_coeff 4 8 lj/cut/coul/long 0.0683 2.8131 # H1-H Ow
pair_coeff 5 5 none # OH-D OH-D
pair_coeff 5 6 none # OH-D HO-D
pair_coeff 5 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # OH-D Hw
pair_coeff 5 8 lj/cut/coul/long/soft 0.1625 3.1427 1.0 # OH-D Ow
pair_coeff 6 6 none # HO-D HO-D
pair_coeff 6 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # HO-D Hw
pair_coeff 6 8 lj/cut/coul/long/soft 0.0000 1.7792 1.0 # HO-D Ow
pair_coeff 7 7 lj/cut/coul/long 0.0000 1.0000 # Hw Hw
pair_coeff 7 8 lj/cut/coul/long 0.0000 1.0000 # Hw Ow
pair_coeff 8 8 lj/cut/coul/long 0.1554 3.1655 # Ow Ow
variable nsteps equal 500000
variable nprint equal ${nsteps}/500
variable ndump equal ${nsteps}/100
# variable nrestart equal ${nsteps}/10
variable temp equal 300.0
variable press equal 1.0
neighbor 2.0 bin
timestep 2.0
# velocity all create ${temp} 12345
thermo_style multi
thermo ${nprint}
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 4 5 a 6
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 500
run 100000
reset_timestep 0
variable dlambda equal -1.0
variable minusdl equal 1.0
variable dq1 equal 0.145*v_dlambda+0.180*v_dlambda
variable dq2 equal 0.180*v_dlambda
variable dq3 equal -0.060*v_dlambda
variable dq4 equal 0.040*v_dlambda-0.060*v_dlambda
variable dq5 equal -0.683*v_dlambda
variable dq6 equal 0.418*v_dlambda
compute cFEP all fep ${temp} &
pair lj/cut/coul/long/soft lambda 2*3 7*8 v_minusdl &
pair lj/cut/coul/long/soft lambda 3 4 v_minusdl &
pair lj/cut/coul/long/soft lambda 4 6 v_dlambda &
pair lj/cut/coul/long/soft lambda 5*6 7*8 v_dlambda &
atom charge 1 v_dq1 &
atom charge 2 v_dq2 &
atom charge 3 v_dq3 &
atom charge 4 v_dq4 &
atom charge 5 v_dq5 &
atom charge 6 v_dq6
fix fFEP all ave/time 1 1 100 c_cFEP[1] c_cFEP[2] file bar10.lmp
# compute cRDF all rdf 100 1 1
# fix fRDF all ave/time 20 100 ${nsteps} c_cRDF file rdf.lammps mode vector
# compute cMSD all msd
# fix fMSD all ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C C H H O H H O
# restart ${nrestart} restart.*.lmp
run ${nsteps}
write_restart restart.*.lmp
write_data data.*.lmp

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,22 @@
# Time-averaged data for fix fFEP
# TimeStep c_cFEP[1] c_cFEP[2]
125000 0.00871763 0.985497
250000 0.00705519 0.98825
375000 0.00568771 0.99052
500000 0.00429336 0.992842
625000 0.00308675 0.994856
750000 0.00169928 0.997181
875000 0.000144366 0.999789
1000000 -0.00112539 1.00193
1125000 -0.0026842 1.00456
1250000 -0.00474318 1.00805
1375000 -0.00571651 1.0097
1500000 -0.0073978 1.01256
1625000 -0.00966334 1.01644
1750000 -0.011584 1.01972
1875000 -0.0148498 1.02535
2000000 -0.0192874 1.03307
2125000 -0.0257123 1.04433
2250000 -0.035944 1.06249
2375000 -0.046465 1.08144
2500000 -0.0558595 1.09863

View File

@ -0,0 +1,142 @@
# created by fftool
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
special_bonds lj/coul 0.0 0.0 0.5
read_data data.0.lmp
# read_restart restart.*.lmp
# reset_timestep 0
pair_style hybrid lj/cut/coul/long 10.0 10.0 &
lj/cut/coul/long/soft 1 0.5 10.0 10.0 10.0
pair_modify tail no
kspace_style pppm 1.0e-4
pair_coeff 1 1 none # C-CO C-CO
pair_coeff 1 2 none # C-CO C-D
pair_coeff 1 3 none # C-CO H-D
pair_coeff 1 4 none # C-CO H-H1
pair_coeff 1 5 none # C-CO OHD-
pair_coeff 1 6 none # C-CO HOD-
pair_coeff 1 7 lj/cut/coul/long 0.0000 1.0000 # C-CO Hw
pair_coeff 1 8 lj/cut/coul/long 0.1013 3.3286 # C-CO Ow
pair_coeff 2 2 none # C-D C-D
pair_coeff 2 3 none # C-D H-D
pair_coeff 2 4 none # C-D H-H1
pair_coeff 2 5 none # C-D OHD-
pair_coeff 2 6 none # C-D HOD-
pair_coeff 2 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # C-D Hw
pair_coeff 2 8 lj/cut/coul/long/soft 0.1013 3.3286 1.0 # C-D Ow
pair_coeff 3 3 none # H-D H-D
pair_coeff 3 4 lj/cut/coul/long/soft 0.0300 2.5000 1.0 # H-D H-H1
pair_coeff 3 5 none # H-D OHD-
pair_coeff 3 6 none # H-D HOD-
pair_coeff 3 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # H-D Hw
pair_coeff 3 8 lj/cut/coul/long/soft 0.0683 2.8131 1.0 # H-D Ow
pair_coeff 4 4 none # H-H1 H-H1
pair_coeff 4 5 none # H-H1 OHD-
pair_coeff 4 6 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # H-H1 HOD-
pair_coeff 4 7 lj/cut/coul/long 0.0000 1.0000 # H-H1 Hw
pair_coeff 4 8 lj/cut/coul/long 0.0683 2.8131 # H-H1 Ow
pair_coeff 5 5 none # OHD- OHD-
pair_coeff 5 6 none # OHD- HOD-
pair_coeff 5 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # OHD- Hw
pair_coeff 5 8 lj/cut/coul/long/soft 0.1625 3.1427 0.0 # OHD- Ow
pair_coeff 6 6 none # HOD- HOD-
pair_coeff 6 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # HOD- Hw
pair_coeff 6 8 lj/cut/coul/long/soft 0.0000 1.7792 0.0 # HOD- Ow
pair_coeff 7 7 lj/cut/coul/long 0.0000 1.0000 # Hw Hw
pair_coeff 7 8 lj/cut/coul/long 0.0000 1.0000 # Hw Ow
pair_coeff 8 8 lj/cut/coul/long 0.1554 3.1655 # Ow Ow
variable nsteps equal 2500000
variable nprint equal ${nsteps}/500
variable ndump equal ${nsteps}/100
# variable nrestart equal ${nsteps}/10
variable temp equal 300.0
variable press equal 1.0
neighbor 2.0 bin
timestep 2.0
velocity all create ${temp} 12345
thermo_style multi
thermo ${nprint}
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 4 5 a 6
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 500
run 250000
reset_timestep 0
variable lambda equal ramp(0.0,1.0)
variable minusl equal 1.0-v_lambda
variable q1 equal 0.145*v_lambda-0.180*(1.0-v_lambda)
variable q2 equal -0.180*(1.0-v_lambda)
variable q3 equal 0.060*(1.0-v_lambda)
variable q4 equal 0.040*v_lambda+0.060*(1.0-v_lambda)
variable q5 equal -0.683*v_lambda
variable q6 equal 0.418*v_lambda
fix fADAPT all adapt/fep 125000 &
pair lj/cut/coul/long/soft lambda 2*3 7*8 v_minusl &
pair lj/cut/coul/long/soft lambda 3 4 v_minusl &
pair lj/cut/coul/long/soft lambda 4 6 v_lambda &
pair lj/cut/coul/long/soft lambda 5*6 7*8 v_lambda &
atom charge 1 v_q1 &
atom charge 2 v_q2 &
atom charge 3 v_q3 &
atom charge 4 v_q4 &
atom charge 5 v_q5 &
atom charge 6 v_q6 &
after yes
variable dlambda equal 0.002
variable minusdl equal -0.002
variable dq1 equal 0.145*v_dlambda+0.180*v_dlambda
variable dq2 equal 0.180*v_dlambda
variable dq3 equal -0.060*v_dlambda
variable dq4 equal 0.040*v_dlambda-0.060*v_dlambda
variable dq5 equal -0.683*v_dlambda
variable dq6 equal 0.418*v_dlambda
compute cFEP all fep ${temp} &
pair lj/cut/coul/long/soft lambda 2*3 7*8 v_minusdl &
pair lj/cut/coul/long/soft lambda 3 4 v_minusdl &
pair lj/cut/coul/long/soft lambda 4 6 v_dlambda &
pair lj/cut/coul/long/soft lambda 5*6 7*8 v_dlambda &
atom charge 1 v_dq1 &
atom charge 2 v_dq2 &
atom charge 3 v_dq3 &
atom charge 4 v_dq4 &
atom charge 5 v_dq5 &
atom charge 6 v_dq6
fix fFEP all ave/time 25 4000 125000 c_cFEP[1] c_cFEP[2] file fdti01.lmp
# compute cRDF all rdf 100 1 1
# fix fRDF all ave/time 20 100 ${nsteps} c_cRDF file rdf.lammps mode vector
compute cMSD all msd
fix fMSD all ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C C H H O H H O
# restart ${nrestart} restart.*.lmp
run ${nsteps}
write_restart restart.*.lmp
write_data data.*.lmp

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,22 @@
# Time-averaged data for fix fFEP
# TimeStep c_cFEP[1] c_cFEP[2]
125000 0.0632017 0.899723
250000 0.0561376 0.910429
375000 0.0471384 0.924277
500000 0.0349153 0.943416
625000 0.0244399 0.960071
750000 0.019451 0.968073
875000 0.0134728 0.977767
1000000 0.011723 0.980628
1125000 0.00945656 0.984351
1250000 0.00790634 0.986905
1375000 0.00564311 0.990642
1500000 0.00376572 0.993753
1625000 0.00238483 0.996053
1750000 0.00085649 0.9986
1875000 -0.000515888 1.0009
2000000 -0.00151569 1.00257
2125000 -0.00325174 1.00549
2250000 -0.00472612 1.00798
2375000 -0.00578231 1.00976
2500000 -0.00729864 1.01233

View File

@ -0,0 +1,142 @@
# created by fftool
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
special_bonds lj/coul 0.0 0.0 0.5
read_data data.1.lmp
# read_restart restart.*.lmp
# reset_timestep 0
pair_style hybrid lj/cut/coul/long 10.0 10.0 &
lj/cut/coul/long/soft 1 0.5 10.0 10.0 10.0
pair_modify tail no
kspace_style pppm 1.0e-4
pair_coeff 1 1 none # CO-C CO-C
pair_coeff 1 2 none # CO-C CD-
pair_coeff 1 3 none # CO-C HD-
pair_coeff 1 4 none # CO-C H1-H
pair_coeff 1 5 none # CO-C OH-D
pair_coeff 1 6 none # CO-C HO-D
pair_coeff 1 7 lj/cut/coul/long 0.0000 1.0000 # CO-C Hw
pair_coeff 1 8 lj/cut/coul/long 0.1013 3.3286 # CO-C Ow
pair_coeff 2 2 none # CD- CD-
pair_coeff 2 3 none # CD- HD-
pair_coeff 2 4 none # CD- H1-H
pair_coeff 2 5 none # CD- OH-D
pair_coeff 2 6 none # CD- HO-D
pair_coeff 2 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # CD- Hw
pair_coeff 2 8 lj/cut/coul/long/soft 0.1013 3.3286 0.0 # CD- Ow
pair_coeff 3 3 none # HD- HD-
pair_coeff 3 4 lj/cut/coul/long/soft 0.0300 2.5000 0.0 # HD- H1-H
pair_coeff 3 5 none # HD- OH-D
pair_coeff 3 6 none # HD- HO-D
pair_coeff 3 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # HD- Hw
pair_coeff 3 8 lj/cut/coul/long/soft 0.0683 2.8131 0.0 # HD- Ow
pair_coeff 4 4 none # H1-H H1-H
pair_coeff 4 5 none # H1-H OH-D
pair_coeff 4 6 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # H1-H HO-D
pair_coeff 4 7 lj/cut/coul/long 0.0000 1.0000 # H1-H Hw
pair_coeff 4 8 lj/cut/coul/long 0.0683 2.8131 # H1-H Ow
pair_coeff 5 5 none # OH-D OH-D
pair_coeff 5 6 none # OH-D HO-D
pair_coeff 5 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # OH-D Hw
pair_coeff 5 8 lj/cut/coul/long/soft 0.1625 3.1427 1.0 # OH-D Ow
pair_coeff 6 6 none # HO-D HO-D
pair_coeff 6 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # HO-D Hw
pair_coeff 6 8 lj/cut/coul/long/soft 0.0000 1.7792 1.0 # HO-D Ow
pair_coeff 7 7 lj/cut/coul/long 0.0000 1.0000 # Hw Hw
pair_coeff 7 8 lj/cut/coul/long 0.0000 1.0000 # Hw Ow
pair_coeff 8 8 lj/cut/coul/long 0.1554 3.1655 # Ow Ow
variable nsteps equal 2500000
variable nprint equal ${nsteps}/500
variable ndump equal ${nsteps}/100
# variable nrestart equal ${nsteps}/10
variable temp equal 300.0
variable press equal 1.0
neighbor 2.0 bin
timestep 2.0
velocity all create ${temp} 12345
thermo_style multi
thermo ${nprint}
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 4 5 a 6
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 500
run 250000
reset_timestep 0
variable lambda equal ramp(1.0,0.0)
variable minusl equal 1.0-v_lambda
variable q1 equal 0.145*v_lambda-0.180*(1.0-v_lambda)
variable q2 equal -0.180*(1.0-v_lambda)
variable q3 equal 0.060*(1.0-v_lambda)
variable q4 equal 0.040*v_lambda+0.060*(1.0-v_lambda)
variable q5 equal -0.683*v_lambda
variable q6 equal 0.418*v_lambda
fix fADAPT all adapt/fep 125000 &
pair lj/cut/coul/long/soft lambda 2*3 7*8 v_minusl &
pair lj/cut/coul/long/soft lambda 3 4 v_minusl &
pair lj/cut/coul/long/soft lambda 4 6 v_lambda &
pair lj/cut/coul/long/soft lambda 5*6 7*8 v_lambda &
atom charge 1 v_q1 &
atom charge 2 v_q2 &
atom charge 3 v_q3 &
atom charge 4 v_q4 &
atom charge 5 v_q5 &
atom charge 6 v_q6 &
after yes
variable dlambda equal -0.002
variable minusdl equal 0.002
variable dq1 equal 0.145*v_dlambda+0.180*v_dlambda
variable dq2 equal 0.180*v_dlambda
variable dq3 equal -0.060*v_dlambda
variable dq4 equal 0.040*v_dlambda-0.060*v_dlambda
variable dq5 equal -0.683*v_dlambda
variable dq6 equal 0.418*v_dlambda
compute cFEP all fep ${temp} &
pair lj/cut/coul/long/soft lambda 2*3 7*8 v_minusdl &
pair lj/cut/coul/long/soft lambda 3 4 v_minusdl &
pair lj/cut/coul/long/soft lambda 4 6 v_dlambda &
pair lj/cut/coul/long/soft lambda 5*6 7*8 v_dlambda &
atom charge 1 v_dq1 &
atom charge 2 v_dq2 &
atom charge 3 v_dq3 &
atom charge 4 v_dq4 &
atom charge 5 v_dq5 &
atom charge 6 v_dq6
fix fFEP all ave/time 25 4000 125000 c_cFEP[1] c_cFEP[2] file fdti10.lmp
# compute cRDF all rdf 100 1 1
# fix fRDF all ave/time 20 100 ${nsteps} c_cRDF file rdf.lammps mode vector
compute cMSD all msd
fix fMSD all ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C C H H O H H O
# restart ${nrestart} restart.*.lmp
run ${nsteps}
write_restart restart.*.lmp
write_data data.*.lmp

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,22 @@
# Time-averaged data for fix fFEP
# TimeStep c_cFEP[1] c_cFEP[2]
125000 0.205889 0.714253
250000 0.164934 0.765732
375000 0.131935 0.809532
500000 0.09771 0.859006
625000 0.0688602 0.903005
750000 0.0356141 0.958919
875000 -0.00180317 1.02349
1000000 -0.0317681 1.08115
1125000 -0.0686221 1.15601
1250000 -0.117495 1.2677
1375000 -0.140627 1.31823
1500000 -0.179791 1.41782
1625000 -0.233651 1.57071
1750000 -0.284345 1.71981
1875000 -0.367527 2.01219
2000000 -0.480496 2.54915
2125000 -0.644218 3.47978
2250000 -0.898797 5.57997
2375000 -1.1516 8.56702
2500000 -1.36364 12.2295

View File

@ -0,0 +1,143 @@
# created by fftool
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
special_bonds lj/coul 0.0 0.0 0.5
read_data data.0.lmp
# read_restart restart.*.lmp
# reset_timestep 0
pair_style hybrid lj/cut/coul/long 10.0 10.0 &
lj/cut/coul/long/soft 1 0.5 10.0 10.0 10.0
pair_modify tail no
kspace_style pppm 1.0e-4
pair_coeff 1 1 none # C-CO C-CO
pair_coeff 1 2 none # C-CO C-D
pair_coeff 1 3 none # C-CO H-D
pair_coeff 1 4 none # C-CO H-H1
pair_coeff 1 5 none # C-CO OHD-
pair_coeff 1 6 none # C-CO HOD-
pair_coeff 1 7 lj/cut/coul/long 0.0000 1.0000 # C-CO Hw
pair_coeff 1 8 lj/cut/coul/long 0.1013 3.3286 # C-CO Ow
pair_coeff 2 2 none # C-D C-D
pair_coeff 2 3 none # C-D H-D
pair_coeff 2 4 none # C-D H-H1
pair_coeff 2 5 none # C-D OHD-
pair_coeff 2 6 none # C-D HOD-
pair_coeff 2 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # C-D Hw
pair_coeff 2 8 lj/cut/coul/long/soft 0.1013 3.3286 1.0 # C-D Ow
pair_coeff 3 3 none # H-D H-D
pair_coeff 3 4 lj/cut/coul/long/soft 0.0300 2.5000 1.0 # H-D H-H1
pair_coeff 3 5 none # H-D OHD-
pair_coeff 3 6 none # H-D HOD-
pair_coeff 3 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # H-D Hw
pair_coeff 3 8 lj/cut/coul/long/soft 0.0683 2.8131 1.0 # H-D Ow
pair_coeff 4 4 none # H-H1 H-H1
pair_coeff 4 5 none # H-H1 OHD-
pair_coeff 4 6 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # H-H1 HOD-
pair_coeff 4 7 lj/cut/coul/long 0.0000 1.0000 # H-H1 Hw
pair_coeff 4 8 lj/cut/coul/long 0.0683 2.8131 # H-H1 Ow
pair_coeff 5 5 none # OHD- OHD-
pair_coeff 5 6 none # OHD- HOD-
pair_coeff 5 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # OHD- Hw
pair_coeff 5 8 lj/cut/coul/long/soft 0.1625 3.1427 0.0 # OHD- Ow
pair_coeff 6 6 none # HOD- HOD-
pair_coeff 6 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # HOD- Hw
pair_coeff 6 8 lj/cut/coul/long/soft 0.0000 1.7792 0.0 # HOD- Ow
pair_coeff 7 7 lj/cut/coul/long 0.0000 1.0000 # Hw Hw
pair_coeff 7 8 lj/cut/coul/long 0.0000 1.0000 # Hw Ow
pair_coeff 8 8 lj/cut/coul/long 0.1554 3.1655 # Ow Ow
variable nsteps equal 25000
variable nprint equal ${nsteps}/500
variable ndump equal ${nsteps}/100
# variable nrestart equal ${nsteps}/10
variable temp equal 300.0
variable press equal 1.0
neighbor 2.0 bin
timestep 2.0
# velocity all create ${temp} 12345
thermo_style multi
thermo ${nprint}
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 4 5 a 6
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 500
run 10000
reset_timestep 0
write_restart init.restart
variable lambda equal ramp(0.0,1.0)
variable minusl equal 1.0-v_lambda
variable q1 equal 0.145*v_lambda-0.180*(1.0-v_lambda)
variable q2 equal -0.180*(1.0-v_lambda)
variable q3 equal 0.060*(1.0-v_lambda)
variable q4 equal 0.040*v_lambda+0.060*(1.0-v_lambda)
variable q5 equal -0.683*v_lambda
variable q6 equal 0.418*v_lambda
fix fADAPT all adapt/fep 125000 &
pair lj/cut/coul/long/soft lambda 2*3 7*8 v_minusl &
pair lj/cut/coul/long/soft lambda 3 4 v_minusl &
pair lj/cut/coul/long/soft lambda 4 6 v_lambda &
pair lj/cut/coul/long/soft lambda 5*6 7*8 v_lambda &
atom charge 1 v_q1 &
atom charge 2 v_q2 &
atom charge 3 v_q3 &
atom charge 4 v_q4 &
atom charge 5 v_q5 &
atom charge 6 v_q6 &
after yes
variable dlambda equal 0.05
variable minusdl equal -0.05
variable dq1 equal 0.145*v_dlambda+0.180*v_dlambda
variable dq2 equal 0.180*v_dlambda
variable dq3 equal -0.060*v_dlambda
variable dq4 equal 0.040*v_dlambda-0.060*v_dlambda
variable dq5 equal -0.683*v_dlambda
variable dq6 equal 0.418*v_dlambda
compute cFEP all fep ${temp} &
pair lj/cut/coul/long/soft lambda 2*3 7*8 v_minusdl &
pair lj/cut/coul/long/soft lambda 3 4 v_minusdl &
pair lj/cut/coul/long/soft lambda 4 6 v_dlambda &
pair lj/cut/coul/long/soft lambda 5*6 7*8 v_dlambda &
atom charge 1 v_dq1 &
atom charge 2 v_dq2 &
atom charge 3 v_dq3 &
atom charge 4 v_dq4 &
atom charge 5 v_dq5 &
atom charge 6 v_dq6
fix fFEP all ave/time 25 4000 125000 c_cFEP[1] c_cFEP[2] file fep01.lmp
# compute cRDF all rdf 100 1 1
# fix fRDF all ave/time 20 100 ${nsteps} c_cRDF file rdf.lammps mode vector
compute cMSD all msd
fix fMSD all ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C C H H O H H O
# restart ${nrestart} restart.*.lmp
run ${nsteps}
write_restart restart.*.lmp
write_data data.*.lmp

View File

@ -0,0 +1,142 @@
# created by fftool
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
special_bonds lj/coul 0.0 0.0 0.5
read_data data.0.lmp
# read_restart restart.*.lmp
# reset_timestep 0
pair_style hybrid lj/cut/coul/long 10.0 10.0 &
lj/cut/coul/long/soft 1 0.5 10.0 10.0 10.0
pair_modify tail no
kspace_style pppm 1.0e-4
pair_coeff 1 1 none # C-CO C-CO
pair_coeff 1 2 none # C-CO C-D
pair_coeff 1 3 none # C-CO H-D
pair_coeff 1 4 none # C-CO H-H1
pair_coeff 1 5 none # C-CO OHD-
pair_coeff 1 6 none # C-CO HOD-
pair_coeff 1 7 lj/cut/coul/long 0.0000 1.0000 # C-CO Hw
pair_coeff 1 8 lj/cut/coul/long 0.1013 3.3286 # C-CO Ow
pair_coeff 2 2 none # C-D C-D
pair_coeff 2 3 none # C-D H-D
pair_coeff 2 4 none # C-D H-H1
pair_coeff 2 5 none # C-D OHD-
pair_coeff 2 6 none # C-D HOD-
pair_coeff 2 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # C-D Hw
pair_coeff 2 8 lj/cut/coul/long/soft 0.1013 3.3286 1.0 # C-D Ow
pair_coeff 3 3 none # H-D H-D
pair_coeff 3 4 lj/cut/coul/long/soft 0.0300 2.5000 1.0 # H-D H-H1
pair_coeff 3 5 none # H-D OHD-
pair_coeff 3 6 none # H-D HOD-
pair_coeff 3 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # H-D Hw
pair_coeff 3 8 lj/cut/coul/long/soft 0.0683 2.8131 1.0 # H-D Ow
pair_coeff 4 4 none # H-H1 H-H1
pair_coeff 4 5 none # H-H1 OHD-
pair_coeff 4 6 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # H-H1 HOD-
pair_coeff 4 7 lj/cut/coul/long 0.0000 1.0000 # H-H1 Hw
pair_coeff 4 8 lj/cut/coul/long 0.0683 2.8131 # H-H1 Ow
pair_coeff 5 5 none # OHD- OHD-
pair_coeff 5 6 none # OHD- HOD-
pair_coeff 5 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # OHD- Hw
pair_coeff 5 8 lj/cut/coul/long/soft 0.1625 3.1427 0.0 # OHD- Ow
pair_coeff 6 6 none # HOD- HOD-
pair_coeff 6 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # HOD- Hw
pair_coeff 6 8 lj/cut/coul/long/soft 0.0000 1.7792 0.0 # HOD- Ow
pair_coeff 7 7 lj/cut/coul/long 0.0000 1.0000 # Hw Hw
pair_coeff 7 8 lj/cut/coul/long 0.0000 1.0000 # Hw Ow
pair_coeff 8 8 lj/cut/coul/long 0.1554 3.1655 # Ow Ow
variable nsteps equal 2500000
variable nprint equal ${nsteps}/500
variable ndump equal ${nsteps}/100
# variable nrestart equal ${nsteps}/10
variable temp equal 300.0
variable press equal 1.0
neighbor 2.0 bin
timestep 2.0
velocity all create ${temp} 12345
thermo_style multi
thermo ${nprint}
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 4 5 a 6
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 500
run 250000
reset_timestep 0
variable lambda equal ramp(0.0,1.0)
variable minusl equal 1.0-v_lambda
variable q1 equal 0.145*v_lambda-0.180*(1.0-v_lambda)
variable q2 equal -0.180*(1.0-v_lambda)
variable q3 equal 0.060*(1.0-v_lambda)
variable q4 equal 0.040*v_lambda+0.060*(1.0-v_lambda)
variable q5 equal -0.683*v_lambda
variable q6 equal 0.418*v_lambda
fix fADAPT all adapt/fep 125000 &
pair lj/cut/coul/long/soft lambda 2*3 7*8 v_minusl &
pair lj/cut/coul/long/soft lambda 3 4 v_minusl &
pair lj/cut/coul/long/soft lambda 4 6 v_lambda &
pair lj/cut/coul/long/soft lambda 5*6 7*8 v_lambda &
atom charge 1 v_q1 &
atom charge 2 v_q2 &
atom charge 3 v_q3 &
atom charge 4 v_q4 &
atom charge 5 v_q5 &
atom charge 6 v_q6 &
after yes
variable dlambda equal 0.05
variable minusdl equal -0.05
variable dq1 equal 0.145*v_dlambda+0.180*v_dlambda
variable dq2 equal 0.180*v_dlambda
variable dq3 equal -0.060*v_dlambda
variable dq4 equal 0.040*v_dlambda-0.060*v_dlambda
variable dq5 equal -0.683*v_dlambda
variable dq6 equal 0.418*v_dlambda
compute cFEP all fep ${temp} &
pair lj/cut/coul/long/soft lambda 2*3 7*8 v_minusdl &
pair lj/cut/coul/long/soft lambda 3 4 v_minusdl &
pair lj/cut/coul/long/soft lambda 4 6 v_dlambda &
pair lj/cut/coul/long/soft lambda 5*6 7*8 v_dlambda &
atom charge 1 v_dq1 &
atom charge 2 v_dq2 &
atom charge 3 v_dq3 &
atom charge 4 v_dq4 &
atom charge 5 v_dq5 &
atom charge 6 v_dq6
fix fFEP all ave/time 25 4000 125000 c_cFEP[1] c_cFEP[2] file fep01.lmp
# compute cRDF all rdf 100 1 1
# fix fRDF all ave/time 20 100 ${nsteps} c_cRDF file rdf.lammps mode vector
compute cMSD all msd
fix fMSD all ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C C H H O H H O
# restart ${nrestart} restart.*.lmp
run ${nsteps}
write_restart restart.*.lmp
write_data data.*.lmp

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,285 @@
LAMMPS (18 Apr 2014-ICMS)
WARNING: OMP_NUM_THREADS environment is not set. (../comm.cpp:100)
using 1 OpenMP thread(s) per MPI task
# created by fftool
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
special_bonds lj/coul 0.0 0.0 0.5
read_data data.0.lmp
orthogonal box = (-13.4438 -13.4438 -13.4438) to (13.4438 13.4438 13.4438)
2 by 2 by 2 MPI processor grid
reading atoms ...
1090 atoms
scanning bonds ...
1 = max bonds/atom
scanning angles ...
10 = max angles/atom
scanning dihedrals ...
12 = max dihedrals/atom
reading bonds ...
729 bonds
reading angles ...
377 angles
reading dihedrals ...
16 dihedrals
5 = max # of 1-2 neighbors
4 = max # of 1-3 neighbors
8 = max # of 1-4 neighbors
9 = max # of special neighbors
# read_restart restart.*.lmp
# reset_timestep 0
pair_style hybrid lj/cut/coul/long 10.0 10.0 lj/cut/coul/long/soft 1 0.5 10.0 10.0 10.0
pair_modify tail no
kspace_style pppm 1.0e-4
pair_coeff 1 1 none # C-CO C-CO
pair_coeff 1 2 none # C-CO C-D
pair_coeff 1 3 none # C-CO H-D
pair_coeff 1 4 none # C-CO H-H1
pair_coeff 1 5 none # C-CO OHD-
pair_coeff 1 6 none # C-CO HOD-
pair_coeff 1 7 lj/cut/coul/long 0.0000 1.0000 # C-CO Hw
pair_coeff 1 8 lj/cut/coul/long 0.1013 3.3286 # C-CO Ow
pair_coeff 2 2 none # C-D C-D
pair_coeff 2 3 none # C-D H-D
pair_coeff 2 4 none # C-D H-H1
pair_coeff 2 5 none # C-D OHD-
pair_coeff 2 6 none # C-D HOD-
pair_coeff 2 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # C-D Hw
pair_coeff 2 8 lj/cut/coul/long/soft 0.1013 3.3286 1.0 # C-D Ow
pair_coeff 3 3 none # H-D H-D
pair_coeff 3 4 lj/cut/coul/long/soft 0.0300 2.5000 1.0 # H-D H-H1
pair_coeff 3 5 none # H-D OHD-
pair_coeff 3 6 none # H-D HOD-
pair_coeff 3 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # H-D Hw
pair_coeff 3 8 lj/cut/coul/long/soft 0.0683 2.8131 1.0 # H-D Ow
pair_coeff 4 4 none # H-H1 H-H1
pair_coeff 4 5 none # H-H1 OHD-
pair_coeff 4 6 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # H-H1 HOD-
pair_coeff 4 7 lj/cut/coul/long 0.0000 1.0000 # H-H1 Hw
pair_coeff 4 8 lj/cut/coul/long 0.0683 2.8131 # H-H1 Ow
pair_coeff 5 5 none # OHD- OHD-
pair_coeff 5 6 none # OHD- HOD-
pair_coeff 5 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # OHD- Hw
pair_coeff 5 8 lj/cut/coul/long/soft 0.1625 3.1427 0.0 # OHD- Ow
pair_coeff 6 6 none # HOD- HOD-
pair_coeff 6 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # HOD- Hw
pair_coeff 6 8 lj/cut/coul/long/soft 0.0000 1.7792 0.0 # HOD- Ow
pair_coeff 7 7 lj/cut/coul/long 0.0000 1.0000 # Hw Hw
pair_coeff 7 8 lj/cut/coul/long 0.0000 1.0000 # Hw Ow
pair_coeff 8 8 lj/cut/coul/long 0.1554 3.1655 # Ow Ow
variable nsteps equal 2500000
variable nprint equal ${nsteps}/500
variable nprint equal 2500000/500
variable ndump equal ${nsteps}/100
variable ndump equal 2500000/100
# variable nrestart equal ${nsteps}/10
variable temp equal 300.0
variable press equal 1.0
neighbor 2.0 bin
timestep 2.0
# velocity all create ${temp} 12345
thermo_style multi
thermo ${nprint}
thermo 5000
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 4 5 a 6
fix fSHAKE all shake 0.0001 20 5000 b 2 4 5 a 6
1 = # of size 2 clusters
0 = # of size 3 clusters
2 = # of size 4 clusters
360 = # of frozen angles
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 500
fix fNPT all npt temp 300 ${temp} 100 iso ${press} ${press} 500
fix fNPT all npt temp 300 300 100 iso ${press} ${press} 500
fix fNPT all npt temp 300 300 100 iso 1 ${press} 500
fix fNPT all npt temp 300 300 100 iso 1 1 500
run 250000
PPPM initialization ...
G vector (1/distance) = 0.270213
grid = 15 15 15
stencil order = 5
estimated absolute RMS force accuracy = 0.0169033
estimated relative force accuracy = 5.09037e-05
using double precision FFTs
3d grid and FFT values/proc = 2744 512
SHAKE stats (type/ave/delta) on step 0
2 1.09 4.99625e-06
4 0.945001 0
5 1 1.58341e-05
6 109.47 0.0012593
Memory usage per processor = 8.38576 Mbytes
---------------- Step 0 ----- CPU = 0.0000 (sec) ----------------
TotEng = -236.7999 KinEng = 0.0000 Temp = 0.0000
PotEng = -236.7999 E_bond = 0.0000 E_angle = 0.3094
E_dihed = 0.0001 E_impro = 0.0000 E_vdwl = -192.9285
E_coul = 19402.7702 E_long = -19446.9511 Press = 406.6448
Volume = 19438.0383
SHAKE stats (type/ave/delta) on step 5000
2 1.09006 8.6791e-07
4 0.945054 0
5 1.00006 3.8639e-06
6 109.47 0.000376049
---------------- Step 5000 ----- CPU = 16.4963 (sec) ----------------
TotEng = -3314.5713 KinEng = 650.7183 Temp = 300.4163
PotEng = -3965.2896 E_bond = 0.0875 E_angle = 1.4729
E_dihed = 0.4749 E_impro = 0.0000 E_vdwl = 774.1508
E_coul = 14908.3382 E_long = -19649.8139 Press = 358.7253
Volume = 10987.1076
SHAKE stats (type/ave/delta) on step 10000
2 1.08998 2.43077e-07
4 0.944985 0
5 0.999985 4.10972e-06
6 109.47 0.000514527
---------------- Step 10000 ----- CPU = 34.0488 (sec) ----------------
TotEng = -3279.4541 KinEng = 689.4215 Temp = 318.2843
PotEng = -3968.8755 E_bond = 0.2040 E_angle = 1.8441
E_dihed = 0.1574 E_impro = 0.0000 E_vdwl = 703.3180
E_coul = 14976.4466 E_long = -19650.8457 Press = -983.5171
Volume = 10927.3524
SHAKE stats (type/ave/delta) on step 15000
2 1.09002 1.0199e-06
4 0.945017 0
5 1.00002 5.08944e-06
6 109.47 0.00042735
---------------- Step 15000 ----- CPU = 53.1541 (sec) ----------------
TotEng = -3368.2566 KinEng = 686.6827 Temp = 317.0199
PotEng = -4054.9393 E_bond = 0.4514 E_angle = 1.5194
E_dihed = 2.4125 E_impro = 0.0000 E_vdwl = 709.7507
E_coul = 14880.7025 E_long = -19649.7759 Press = -1277.1945
Volume = 10817.2760
SHAKE stats (type/ave/delta) on step 20000
2 1.09005 3.84509e-06
4 0.945044 0
5 1.00005 1.15463e-05
6 109.47 0.000893174
---------------- Step 20000 ----- CPU = 67.0658 (sec) ----------------
TotEng = -3401.0672 KinEng = 639.1818 Temp = 295.0902
PotEng = -4040.2490 E_bond = 0.8254 E_angle = 4.2422
E_dihed = 0.6223 E_impro = 0.0000 E_vdwl = 777.7799
E_coul = 14826.3237 E_long = -19650.0425 Press = -268.5153
Volume = 10911.3140
SHAKE stats (type/ave/delta) on step 25000
2 1.09 3.16085e-06
4 0.945002 0
5 1 6.82105e-06
6 109.47 0.000683982
---------------- Step 25000 ----- CPU = 83.1242 (sec) ----------------
TotEng = -3331.9648 KinEng = 658.5638 Temp = 304.0383
PotEng = -3990.5287 E_bond = 0.0504 E_angle = 4.8757
E_dihed = 0.9548 E_impro = 0.0000 E_vdwl = 742.1312
E_coul = 14907.2989 E_long = -19645.8396 Press = -803.8894
Volume = 11243.5807
SHAKE stats (type/ave/delta) on step 30000
2 1.08995 5.90921e-08
4 0.944954 0
5 0.999952 5.01043e-06
6 109.47 0.000369665
---------------- Step 30000 ----- CPU = 96.1911 (sec) ----------------
TotEng = -3321.5197 KinEng = 681.8132 Temp = 314.7718
PotEng = -4003.3328 E_bond = 0.5162 E_angle = 6.2281
E_dihed = 1.2145 E_impro = 0.0000 E_vdwl = 767.4719
E_coul = 14871.6899 E_long = -19650.4535 Press = 63.7724
Volume = 10913.5775
SHAKE stats (type/ave/delta) on step 35000
2 1.08996 3.93868e-07
4 0.944963 0
5 0.999961 4.01003e-06
6 109.47 0.000492251
---------------- Step 35000 ----- CPU = 111.1547 (sec) ----------------
TotEng = -3376.3100 KinEng = 652.5134 Temp = 301.2450
PotEng = -4028.8234 E_bond = 0.6940 E_angle = 3.8874
E_dihed = 0.6796 E_impro = 0.0000 E_vdwl = 811.7936
E_coul = 14803.9468 E_long = -19649.8247 Press = 957.3673
Volume = 10810.8812
SHAKE stats (type/ave/delta) on step 40000
2 1.08993 1.84911e-06
4 0.944937 0
5 0.999934 4.48159e-06
6 109.47 0.000423379
---------------- Step 40000 ----- CPU = 127.8336 (sec) ----------------
TotEng = -3342.5048 KinEng = 641.0827 Temp = 295.9678
PotEng = -3983.5875 E_bond = 1.4592 E_angle = 3.3579
E_dihed = 0.1250 E_impro = 0.0000 E_vdwl = 714.5090
E_coul = 14945.5914 E_long = -19648.6299 Press = -1146.3509
Volume = 10993.6945
SHAKE stats (type/ave/delta) on step 45000
2 1.09004 2.52413e-06
4 0.945038 0
5 1.00004 7.30752e-06
6 109.47 0.000660905
---------------- Step 45000 ----- CPU = 141.7030 (sec) ----------------
TotEng = -3368.5223 KinEng = 640.7082 Temp = 295.7949
PotEng = -4009.2305 E_bond = 0.0226 E_angle = 2.7492
E_dihed = 0.8929 E_impro = 0.0000 E_vdwl = 820.3025
E_coul = 14817.6347 E_long = -19650.8325 Press = 388.8981
Volume = 11292.9397
SHAKE stats (type/ave/delta) on step 50000
2 1.09006 2.89749e-06
4 0.945051 0
5 1.00006 6.08927e-06
6 109.47 0.000779845
---------------- Step 50000 ----- CPU = 156.2936 (sec) ----------------
TotEng = -3253.6028 KinEng = 683.0893 Temp = 315.3610
PotEng = -3936.6921 E_bond = 0.4303 E_angle = 5.1773
E_dihed = 1.1538 E_impro = 0.0000 E_vdwl = 722.3282
E_coul = 14982.9162 E_long = -19648.6979 Press = -338.7238
Volume = 11230.1507
SHAKE stats (type/ave/delta) on step 55000
2 1.09007 9.95472e-07
4 0.945062 0
5 1.00007 3.93786e-06
6 109.47 0.000436645
---------------- Step 55000 ----- CPU = 176.2674 (sec) ----------------
TotEng = -3332.0820 KinEng = 662.9270 Temp = 306.0527
PotEng = -3995.0090 E_bond = 0.8395 E_angle = 1.9527
E_dihed = 1.0770 E_impro = 0.0000 E_vdwl = 750.6359
E_coul = 14899.3802 E_long = -19648.8943 Press = -498.7755
Volume = 11082.6763
SHAKE stats (type/ave/delta) on step 60000
2 1.09006 7.79087e-08
4 0.945052 1.22125e-15
5 1.00006 3.87615e-06
6 109.47 0.000355506
---------------- Step 60000 ----- CPU = 193.3379 (sec) ----------------
TotEng = -3296.5834 KinEng = 679.6746 Temp = 313.7845
PotEng = -3976.2580 E_bond = 2.4610 E_angle = 2.9231
E_dihed = 0.8719 E_impro = 0.0000 E_vdwl = 803.1348
E_coul = 14863.3163 E_long = -19648.9651 Press = 1009.2429
Volume = 11125.1736
SHAKE stats (type/ave/delta) on step 65000
2 1.09003 1.05995e-06
4 0.945023 1.22125e-15
5 1.00003 4.3471e-06
6 109.47 0.000511365
---------------- Step 65000 ----- CPU = 210.3194 (sec) ----------------
TotEng = -3322.0496 KinEng = 666.5091 Temp = 307.7064
PotEng = -3988.5587 E_bond = 0.8728 E_angle = 5.7476
E_dihed = 0.9783 E_impro = 0.0000 E_vdwl = 743.8165
E_coul = 14910.0426 E_long = -19650.0166 Press = -421.7167
Volume = 10902.5998
SHAKE stats (type/ave/delta) on step 70000
2 1.0901 1.09831e-06
4 0.94509 0
5 1.0001 3.59309e-06
6 109.47 0.000361318
---------------- Step 70000 ----- CPU = 228.7878 (sec) ----------------
TotEng = -3395.7912 KinEng = 667.4514 Temp = 308.1414
PotEng = -4063.2426 E_bond = 0.6965 E_angle = 5.1350
E_dihed

View File

@ -0,0 +1,287 @@
LAMMPS (18 Apr 2014-ICMS)
WARNING: OMP_NUM_THREADS environment is not set. (../comm.cpp:100)
using 1 OpenMP thread(s) per MPI task
package omp *
using multi-threaded neighbor list subroutines
prefer double precision OpenMP force kernels
# created by fftool
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
special_bonds lj/coul 0.0 0.0 0.5
read_data data.0.lmp
orthogonal box = (-13.4438 -13.4438 -13.4438) to (13.4438 13.4438 13.4438)
2 by 2 by 2 MPI processor grid
reading atoms ...
1090 atoms
scanning bonds ...
1 = max bonds/atom
scanning angles ...
10 = max angles/atom
scanning dihedrals ...
12 = max dihedrals/atom
reading bonds ...
729 bonds
reading angles ...
377 angles
reading dihedrals ...
16 dihedrals
5 = max # of 1-2 neighbors
4 = max # of 1-3 neighbors
8 = max # of 1-4 neighbors
9 = max # of special neighbors
# read_restart restart.*.lmp
# reset_timestep 0
pair_style hybrid lj/cut/coul/long 10.0 10.0 lj/cut/coul/long/soft 1 0.5 10.0 10.0 10.0
pair_modify tail no
kspace_style pppm 1.0e-4
pair_coeff 1 1 none # C-CO C-CO
pair_coeff 1 2 none # C-CO C-D
pair_coeff 1 3 none # C-CO H-D
pair_coeff 1 4 none # C-CO H-H1
pair_coeff 1 5 none # C-CO OHD-
pair_coeff 1 6 none # C-CO HOD-
pair_coeff 1 7 lj/cut/coul/long 0.0000 1.0000 # C-CO Hw
pair_coeff 1 8 lj/cut/coul/long 0.1013 3.3286 # C-CO Ow
pair_coeff 2 2 none # C-D C-D
pair_coeff 2 3 none # C-D H-D
pair_coeff 2 4 none # C-D H-H1
pair_coeff 2 5 none # C-D OHD-
pair_coeff 2 6 none # C-D HOD-
pair_coeff 2 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # C-D Hw
pair_coeff 2 8 lj/cut/coul/long/soft 0.1013 3.3286 1.0 # C-D Ow
pair_coeff 3 3 none # H-D H-D
pair_coeff 3 4 lj/cut/coul/long/soft 0.0300 2.5000 1.0 # H-D H-H1
pair_coeff 3 5 none # H-D OHD-
pair_coeff 3 6 none # H-D HOD-
pair_coeff 3 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # H-D Hw
pair_coeff 3 8 lj/cut/coul/long/soft 0.0683 2.8131 1.0 # H-D Ow
pair_coeff 4 4 none # H-H1 H-H1
pair_coeff 4 5 none # H-H1 OHD-
pair_coeff 4 6 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # H-H1 HOD-
pair_coeff 4 7 lj/cut/coul/long 0.0000 1.0000 # H-H1 Hw
pair_coeff 4 8 lj/cut/coul/long 0.0683 2.8131 # H-H1 Ow
pair_coeff 5 5 none # OHD- OHD-
pair_coeff 5 6 none # OHD- HOD-
pair_coeff 5 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # OHD- Hw
pair_coeff 5 8 lj/cut/coul/long/soft 0.1625 3.1427 0.0 # OHD- Ow
pair_coeff 6 6 none # HOD- HOD-
pair_coeff 6 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # HOD- Hw
pair_coeff 6 8 lj/cut/coul/long/soft 0.0000 1.7792 0.0 # HOD- Ow
pair_coeff 7 7 lj/cut/coul/long 0.0000 1.0000 # Hw Hw
pair_coeff 7 8 lj/cut/coul/long 0.0000 1.0000 # Hw Ow
pair_coeff 8 8 lj/cut/coul/long 0.1554 3.1655 # Ow Ow
variable nsteps equal 2500000
variable nprint equal ${nsteps}/500
variable nprint equal 2500000/500
variable ndump equal ${nsteps}/100
variable ndump equal 2500000/100
# variable nrestart equal ${nsteps}/10
variable temp equal 300.0
variable press equal 1.0
neighbor 2.0 bin
timestep 2.0
# velocity all create ${temp} 12345
thermo_style multi
thermo ${nprint}
thermo 5000
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 4 5 a 6
fix fSHAKE all shake 0.0001 20 5000 b 2 4 5 a 6
1 = # of size 2 clusters
0 = # of size 3 clusters
2 = # of size 4 clusters
360 = # of frozen angles
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 500
fix fNPT all npt temp 300 ${temp} 100 iso ${press} ${press} 500
fix fNPT all npt temp 300 300 100 iso ${press} ${press} 500
fix fNPT all npt temp 300 300 100 iso 1 ${press} 500
fix fNPT all npt temp 300 300 100 iso 1 1 500
run 250000
PPPM initialization ...
G vector (1/distance) = 0.270213
grid = 15 15 15
stencil order = 5
estimated absolute RMS force accuracy = 0.0169033
estimated relative force accuracy = 5.09037e-05
using double precision FFTs
3d grid and FFT values/proc = 2744 512
Hybrid pair style last /omp style lj/cut/coul/long/soft
Last active /omp style is kspace_style pppm/omp
SHAKE stats (type/ave/delta) on step 0
2 1.09 4.99625e-06
4 0.945001 0
5 1 1.58341e-05
6 109.47 0.0012593
Memory usage per processor = 8.38667 Mbytes
---------------- Step 0 ----- CPU = 0.0000 (sec) ----------------
TotEng = -236.7999 KinEng = 0.0000 Temp = 0.0000
PotEng = -236.7999 E_bond = 0.0000 E_angle = 0.3094
E_dihed = 0.0001 E_impro = 0.0000 E_vdwl = -192.9285
E_coul = 19402.7702 E_long = -19446.9511 Press = 406.6448
Volume = 19438.0383
SHAKE stats (type/ave/delta) on step 5000
2 1.09009 2.87398e-07
4 0.945079 0
5 1.00008 3.9811e-06
6 109.47 0.000369321
---------------- Step 5000 ----- CPU = 11.8253 (sec) ----------------
TotEng = -3307.2904 KinEng = 658.6061 Temp = 304.0578
PotEng = -3965.8965 E_bond = 0.0083 E_angle = 0.6669
E_dihed = 0.2569 E_impro = 0.0000 E_vdwl = 734.5556
E_coul = 14945.5504 E_long = -19646.9346 Press = -401.3652
Volume = 10990.0279
SHAKE stats (type/ave/delta) on step 10000
2 1.08991 8.10164e-07
4 0.944925 0
5 0.999922 7.76553e-06
6 109.47 0.000694232
---------------- Step 10000 ----- CPU = 23.7032 (sec) ----------------
TotEng = -3381.5340 KinEng = 646.2577 Temp = 298.3570
PotEng = -4027.7917 E_bond = 0.2113 E_angle = 1.2103
E_dihed = 0.8831 E_impro = 0.0000 E_vdwl = 807.4931
E_coul = 14812.7871 E_long = -19650.3766 Press = 1168.7845
Volume = 10591.5096
SHAKE stats (type/ave/delta) on step 15000
2 1.08988 6.96695e-07
4 0.944894 0
5 0.999889 7.87681e-06
6 109.47 0.000849755
---------------- Step 15000 ----- CPU = 37.0159 (sec) ----------------
TotEng = -3391.4919 KinEng = 640.6926 Temp = 295.7878
PotEng = -4032.1845 E_bond = 0.3271 E_angle = 2.4027
E_dihed = 0.2102 E_impro = 0.0000 E_vdwl = 833.3641
E_coul = 14781.8806 E_long = -19650.3692 Press = 1246.1602
Volume = 10806.8631
SHAKE stats (type/ave/delta) on step 20000
2 1.09009 5.35061e-07
4 0.945076 0
5 1.00008 5.8947e-06
6 109.47 0.000481652
---------------- Step 20000 ----- CPU = 49.5427 (sec) ----------------
TotEng = -3340.9663 KinEng = 649.0129 Temp = 299.6290
PotEng = -3989.9791 E_bond = 0.2338 E_angle = 1.4422
E_dihed = 0.7066 E_impro = 0.0000 E_vdwl = 823.4408
E_coul = 14835.4457 E_long = -19651.2482 Press = 1369.0812
Volume = 10984.3808
SHAKE stats (type/ave/delta) on step 25000
2 1.09003 1.29003e-06
4 0.945026 0
5 1.00003 5.08541e-06
6 109.47 0.000543653
---------------- Step 25000 ----- CPU = 61.5352 (sec) ----------------
TotEng = -3350.1301 KinEng = 628.6694 Temp = 290.2370
PotEng = -3978.7994 E_bond = 0.5250 E_angle = 3.3362
E_dihed = 1.1997 E_impro = 0.0000 E_vdwl = 753.6460
E_coul = 14908.8825 E_long = -19646.3890 Press = -107.3383
Volume = 10858.4423
SHAKE stats (type/ave/delta) on step 30000
2 1.09012 2.15723e-06
4 0.945104 0
5 1.00011 7.30223e-06
6 109.47 0.000722195
---------------- Step 30000 ----- CPU = 75.1590 (sec) ----------------
TotEng = -3291.1533 KinEng = 662.9659 Temp = 306.0706
PotEng = -3954.1192 E_bond = 0.0386 E_angle = 3.4541
E_dihed = 0.3158 E_impro = 0.0000 E_vdwl = 745.0479
E_coul = 14947.2099 E_long = -19650.1856 Press = -296.2930
Volume = 11301.5139
SHAKE stats (type/ave/delta) on step 35000
2 1.08998 1.28162e-06
4 0.944981 0
5 0.999982 8.25561e-06
6 109.47 0.000790911
---------------- Step 35000 ----- CPU = 86.8870 (sec) ----------------
TotEng = -3370.6689 KinEng = 637.0075 Temp = 294.0864
PotEng = -4007.6763 E_bond = 0.2585 E_angle = 2.7426
E_dihed = 0.2288 E_impro = 0.0000 E_vdwl = 766.3429
E_coul = 14871.1344 E_long = -19648.3835 Press = -434.1739
Volume = 11077.2342
SHAKE stats (type/ave/delta) on step 40000
2 1.09005 3.30799e-07
4 0.945044 0
5 1.00005 4.39187e-06
6 109.47 0.000367325
---------------- Step 40000 ----- CPU = 102.1257 (sec) ----------------
TotEng = -3318.0442 KinEng = 694.3948 Temp = 320.5803
PotEng = -4012.4390 E_bond = 0.3986 E_angle = 4.3594
E_dihed = 1.0551 E_impro = 0.0000 E_vdwl = 716.6662
E_coul = 14914.0530 E_long = -19648.9712 Press = -1177.4188
Volume = 11087.6184
SHAKE stats (type/ave/delta) on step 45000
2 1.08997 7.63663e-06
4 0.944973 0
5 0.999973 2.30441e-05
6 109.47 0.00155851
---------------- Step 45000 ----- CPU = 119.3648 (sec) ----------------
TotEng = -3371.5332 KinEng = 660.1645 Temp = 304.7773
PotEng = -4031.6977 E_bond = 1.3352 E_angle = 2.6760
E_dihed = 1.0255 E_impro = 0.0000 E_vdwl = 734.4320
E_coul = 14878.0777 E_long = -19649.2440 Press = -978.3955
Volume = 11108.5059
SHAKE stats (type/ave/delta) on step 50000
2 1.09015 1.28872e-06
4 0.945128 0
5 1.00014 5.27234e-06
6 109.47 0.000408939
---------------- Step 50000 ----- CPU = 132.4029 (sec) ----------------
TotEng = -3345.1117 KinEng = 660.2366 Temp = 304.8106
PotEng = -4005.3482 E_bond = 0.9878 E_angle = 3.3427
E_dihed = 0.4740 E_impro = 0.0000 E_vdwl = 763.3350
E_coul = 14874.5886 E_long = -19648.0763 Press = 257.8221
Volume = 10750.8743
SHAKE stats (type/ave/delta) on step 55000
2 1.09006 1.79853e-06
4 0.945055 0
5 1.00006 8.96473e-06
6 109.47 0.000604022
---------------- Step 55000 ----- CPU = 144.3965 (sec) ----------------
TotEng = -3329.0112 KinEng = 639.1918 Temp = 295.0949
PotEng = -3968.2030 E_bond = 0.3779 E_angle = 4.2627
E_dihed = 0.1647 E_impro = 0.0000 E_vdwl = 775.0143
E_coul = 14900.7482 E_long = -19648.7708 Press = 384.3466
Volume = 10943.9465
SHAKE stats (type/ave/delta) on step 60000
2 1.09006 5.76725e-07
4 0.94505 0
5 1.00005 4.15938e-06
6 109.47 0.000456895
---------------- Step 60000 ----- CPU = 158.9285 (sec) ----------------
TotEng = -3337.6373 KinEng = 686.3391 Temp = 316.8613
PotEng = -4023.9764 E_bond = 0.1951 E_angle = 3.2295
E_dihed = 0.2227 E_impro = 0.0000 E_vdwl = 803.8460
E_coul = 14822.5916 E_long = -19654.0613 Press = 343.2839
Volume = 11053.5476
SHAKE stats (type/ave/delta) on step 65000
2 1.08992 1.92627e-06
4 0.944928 0
5 0.999925 9.65284e-06
6 109.47 0.000982037
---------------- Step 65000 ----- CPU = 171.6679 (sec) ----------------
TotEng = -3286.1431 KinEng = 655.6296 Temp = 302.6837
PotEng = -3941.7726 E_bond = 0.8912 E_angle = 5.2439
E_dihed = 1.0392 E_impro = 0.0000 E_vdwl = 723.8338
E_coul = 14977.3056 E_long = -19650.0863 Press = -610.8426
Volume = 11125.9793
SHAKE stats (type/ave/delta) on step 70000
2 1.08995 3.49023e-07
4 0.944954 0
5 0.999952 3.62092e-06
6 109.47 0.000402967
---------------- Step 70000 ----- CPU

View File

@ -0,0 +1,22 @@
# Time-averaged data for fix fFEP
# TimeStep c_cFEP[1] c_cFEP[2]
125000 1.6302 0.0823708
250000 1.42237 0.11283
375000 1.1798 0.168088
500000 0.869449 0.279336
625000 0.610803 0.412455
750000 0.489597 0.486584
875000 0.34439 0.601048
1000000 0.303411 0.640696
1125000 0.244627 0.70206
1250000 0.204737 0.746389
1375000 0.144716 0.816302
1500000 0.0941156 0.880968
1625000 0.0575351 0.933564
1750000 0.0171831 0.992497
1875000 -0.018964 1.05078
2000000 -0.0457397 1.09726
2125000 -0.0904898 1.17888
2250000 -0.128205 1.25226
2375000 -0.155891 1.31122
2500000 -0.194237 1.3965

View File

@ -0,0 +1,142 @@
# created by fftool
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
special_bonds lj/coul 0.0 0.0 0.5
read_data data.1.lmp
# read_restart restart.*.lmp
# reset_timestep 0
pair_style hybrid lj/cut/coul/long 10.0 10.0 &
lj/cut/coul/long/soft 1 0.5 10.0 10.0 10.0
pair_modify tail no
kspace_style pppm 1.0e-4
pair_coeff 1 1 none # CO-C CO-C
pair_coeff 1 2 none # CO-C CD-
pair_coeff 1 3 none # CO-C HD-
pair_coeff 1 4 none # CO-C H1-H
pair_coeff 1 5 none # CO-C OH-D
pair_coeff 1 6 none # CO-C HO-D
pair_coeff 1 7 lj/cut/coul/long 0.0000 1.0000 # CO-C Hw
pair_coeff 1 8 lj/cut/coul/long 0.1013 3.3286 # CO-C Ow
pair_coeff 2 2 none # CD- CD-
pair_coeff 2 3 none # CD- HD-
pair_coeff 2 4 none # CD- H1-H
pair_coeff 2 5 none # CD- OH-D
pair_coeff 2 6 none # CD- HO-D
pair_coeff 2 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # CD- Hw
pair_coeff 2 8 lj/cut/coul/long/soft 0.1013 3.3286 0.0 # CD- Ow
pair_coeff 3 3 none # HD- HD-
pair_coeff 3 4 lj/cut/coul/long/soft 0.0300 2.5000 0.0 # HD- H1-H
pair_coeff 3 5 none # HD- OH-D
pair_coeff 3 6 none # HD- HO-D
pair_coeff 3 7 lj/cut/coul/long/soft 0.0000 1.0000 0.0 # HD- Hw
pair_coeff 3 8 lj/cut/coul/long/soft 0.0683 2.8131 0.0 # HD- Ow
pair_coeff 4 4 none # H1-H H1-H
pair_coeff 4 5 none # H1-H OH-D
pair_coeff 4 6 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # H1-H HO-D
pair_coeff 4 7 lj/cut/coul/long 0.0000 1.0000 # H1-H Hw
pair_coeff 4 8 lj/cut/coul/long 0.0683 2.8131 # H1-H Ow
pair_coeff 5 5 none # OH-D OH-D
pair_coeff 5 6 none # OH-D HO-D
pair_coeff 5 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # OH-D Hw
pair_coeff 5 8 lj/cut/coul/long/soft 0.1625 3.1427 1.0 # OH-D Ow
pair_coeff 6 6 none # HO-D HO-D
pair_coeff 6 7 lj/cut/coul/long/soft 0.0000 1.0000 1.0 # HO-D Hw
pair_coeff 6 8 lj/cut/coul/long/soft 0.0000 1.7792 1.0 # HO-D Ow
pair_coeff 7 7 lj/cut/coul/long 0.0000 1.0000 # Hw Hw
pair_coeff 7 8 lj/cut/coul/long 0.0000 1.0000 # Hw Ow
pair_coeff 8 8 lj/cut/coul/long 0.1554 3.1655 # Ow Ow
variable nsteps equal 2500000
variable nprint equal ${nsteps}/500
variable ndump equal ${nsteps}/100
# variable nrestart equal ${nsteps}/10
variable temp equal 300.0
variable press equal 1.0
neighbor 2.0 bin
timestep 2.0
velocity all create ${temp} 12345
thermo_style multi
thermo ${nprint}
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 4 5 a 6
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 500
run 250000
reset_timestep 0
variable lambda equal ramp(1.0,0.0)
variable minusl equal 1.0-v_lambda
variable q1 equal 0.145*v_lambda-0.180*(1.0-v_lambda)
variable q2 equal -0.180*(1.0-v_lambda)
variable q3 equal 0.060*(1.0-v_lambda)
variable q4 equal 0.040*v_lambda+0.060*(1.0-v_lambda)
variable q5 equal -0.683*v_lambda
variable q6 equal 0.418*v_lambda
fix fADAPT all adapt/fep 125000 &
pair lj/cut/coul/long/soft lambda 2*3 7*8 v_minusl &
pair lj/cut/coul/long/soft lambda 3 4 v_minusl &
pair lj/cut/coul/long/soft lambda 4 6 v_lambda &
pair lj/cut/coul/long/soft lambda 5*6 7*8 v_lambda &
atom charge 1 v_q1 &
atom charge 2 v_q2 &
atom charge 3 v_q3 &
atom charge 4 v_q4 &
atom charge 5 v_q5 &
atom charge 6 v_q6 &
after yes
variable dlambda equal -0.05
variable minusdl equal 0.05
variable dq1 equal 0.145*v_dlambda+0.180*v_dlambda
variable dq2 equal 0.180*v_dlambda
variable dq3 equal -0.060*v_dlambda
variable dq4 equal 0.040*v_dlambda-0.060*v_dlambda
variable dq5 equal -0.683*v_dlambda
variable dq6 equal 0.418*v_dlambda
compute cFEP all fep ${temp} &
pair lj/cut/coul/long/soft lambda 2*3 7*8 v_minusdl &
pair lj/cut/coul/long/soft lambda 3 4 v_minusdl &
pair lj/cut/coul/long/soft lambda 4 6 v_dlambda &
pair lj/cut/coul/long/soft lambda 5*6 7*8 v_dlambda &
atom charge 1 v_dq1 &
atom charge 2 v_dq2 &
atom charge 3 v_dq3 &
atom charge 4 v_dq4 &
atom charge 5 v_dq5 &
atom charge 6 v_dq6
fix fFEP all ave/time 25 4000 125000 c_cFEP[1] c_cFEP[2] file fep10.lmp
# compute cRDF all rdf 100 1 1
# fix fRDF all ave/time 20 100 ${nsteps} c_cRDF file rdf.lammps mode vector
compute cMSD all msd
fix fMSD all ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C C H H O H H O
# restart ${nrestart} restart.*.lmp
run ${nsteps}
write_restart restart.*.lmp
write_data data.*.lmp

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,49 @@
# CC-CO.ff (kJ/mol, A, deg)
ATOMS
# type m/uma q/e pot pars
# ethane OPLS-AA
C-CO CT 12.011 -0.180 lj 3.50 0.27614
C-D CT 12.011 -0.180 lj 3.50 0.27614
CD- CT 12.011 0.000 lj 3.50 0.27614
H-D HC 1.008 0.060 lj 2.50 0.12552
HD- HC 1.008 0.000 lj 2.50 0.12552
H-H1 HC 1.008 0.060 lj 2.50 0.12552
# methanol OPLS-AA
CO-C CT 12.011 0.145 lj 3.50 0.27614
H1-H HC 1.008 0.040 lj 2.50 0.12552
OH-D OH 15.999 -0.683 lj 3.12 0.71128
OHD- OH 15.999 0.000 lj 3.12 0.71128
HO-D HO 1.008 0.418 lj 1.00 0.00000
HOD- HO 1.008 0.000 lj 1.00 0.00000
BONDS
# i j pot re/A ka/kJmol-1
# alkanes OPLS-AA JACS 118(1996)11225, JPC 100(1996)18010
HC CT cons 1.090 2845.12
CT CT harm 1.529 2242.62
# alcohols OPLS-AA JACS 118(1996)11225, JPC 100(1996)18010
CT OH harm 1.410 2677.76
HO OH cons 0.945 4627.50
ANGLES
# i j k pot th/deg ka/kjmol-1
# alkanes OPLS-AA JACS 118(1996)11225, JPC 100(1996)18010
HC CT HC harm 107.80 276.14
HC CT CT harm 110.70 313.80
CT CT CT harm 112.70 488.27
# alcohols OPLS-AA JACS 118(1996)11225, JPC 100(1996)18010
CT CT OH harm 0.0 418.40
HC CT OH harm 109.50 292.88
CT OH HO harm 108.50 460.24
DIHEDRALS
# i j k l pot v1 v2 v3 v4
# alkanes OPLS-AA JACS 118: 11225 (1996); J Phys Chem 100: 18010 (1996)
HC CT CT HC opls 0.0000 0.0000 1.3305 0.0000
HC CT CT CT opls 0.0000 0.0000 1.5313 0.0000
CT CT CT CT opls 7.2800 -0.6569 1.1673 0.0000
# alcohols OPLS-AA JACS 118(1996)11225, AMBER98 (OCCO) 117(1995)5179
HC CT OH HO opls 0.0000 0.0000 1.8828 0.0000
CT CT OH HO opls 0.0000 0.0000 0.0000 0.0000
HC CT CT OH opls 0.0000 0.0000 0.0000 0.0000

View File

@ -0,0 +1,14 @@
CC->CO
1 C-CO
2 C-D 1 1.529
3 H-D 2 1.090 1 109.5
4 H-D 2 1.090 1 109.5 3 120.0
5 H-D 2 1.090 1 109.5 3 240.0
6 H-H1 1 1.090 2 109.5 3 180.0
7 H-H1 1 1.090 2 109.5 4 180.0
8 H-H1 1 1.090 2 109.5 5 180.0
9 OHD- 1 1.410 2 1.0 3 0.0
10 HOD- 9 0.945 1 109.5 2 180.0
CC-CO.ff

View File

@ -0,0 +1,14 @@
CO->CC
1 CO-C
2 CD- 1 1.529
3 HD- 2 1.090 1 109.5
4 HD- 2 1.090 1 109.5 3 120.0
5 HD- 2 1.090 1 109.5 3 240.0
6 H1-H 1 1.090 2 109.5 3 180.0
7 H1-H 1 1.090 2 109.5 4 180.0
8 H1-H 1 1.090 2 109.5 5 180.0
9 OH-D 1 1.410 2 1.0 3 0.0
10 HO-D 9 0.945 1 109.5 2 180.0
CC-CO.ff

View File

@ -0,0 +1,15 @@
Ethane to Methanol Molecule Files
=================================
Files necessary to build starting configurations using the
[fftool.py](http://www.github.com/agiliopadua/fftool) utility and
[Packmol](http://www.ime.unicamp.br/~martinez/packmol/), for example:
fftool.py 1 CC-CO.zmat 360 spce.zmat -r 50
packmol < pack.in
fftool.py 1 CC-CO.zmat 360 spce.zmat -r 50 -l -a
The resulting `in.lmp` files have to be modified to introduce the free
energy route.

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,14 @@
# SPC/E water (kJ/mol, A, deg)
ATOMS
# type m/uma q/e pot pars
Ow Ow 15.999 -0.8476 lj 3.1655 0.6503
Hw Hw 1.008 0.4238 lj 1.0000 0.0000
BONDS
# i j pot re/A ka/kJmol-1
Ow Hw cons 1.000 4331.53
ANGLES
# i j k pot th/deg ka/kjmol-1
Hw Ow Hw cons 109.47 317.57

View File

@ -0,0 +1,7 @@
SPCE
Hw
Ow 1 1.000
Hw 2 1.000 1 109.47
spce.ff

View File

@ -0,0 +1,57 @@
Methane to Tetrafluoromethane in Water
======================================
Example calculation of the difference in free energy of hydration upon
transforming methane into tetrafluoromethane with LAMMPS using
*compute fep* and *fix adapt*.
Methane and tetrafluoromethane are represented by the OPLS-AA force
field (1 molecule). Water is represented by the 3-site SPC/E model
(360 molecules).
The strategy used to perform the alchemical transformation is the
following:
* The dual topology approach is used, therefore all the atoms of
methane and perfluorommethane are present throughout the simulation,
only some of them are dummy sites at the endpoints of the
transformation. Masses and intramolecular terms (bond lengths,
angles, dihedrals) are not changed.
* Interactions of sites that are being created (from dummy sites) or
deleted (to become dummy sites) are treated using soft-core verions
of the Lennard-Jones and Coulomb potentials (*pair
lj/cut/coul/long/soft*) in order to avoid singularities. The
exponent of the coupling parameter lambda in the soft-core pair
potentials was in this example n = 1.
* Eletrostatic charges that are modified are varied linearly from the
initial to the final values. This keeps the overall charge of the
molecule constant, which is good for the long range electrostatics
(the coupling parameter lambda has no effect on the kspace terms).
The following directories contain input files and results for
calculations using Bennet's acceptance ratio (BAR) method:
* `bar01` -- Calculation using BAR, 1-step transformation of a CH4
molecule into CF4. Results in `bar01.lmp`
* `bar10` -- Calculation using BAR, 1-step transformation of a
CF4 molecule into CH4. Results in `bar10.lmp`
The Python script `bar.py` found in the `tools` directory can
be used to calculate the free-energy difference corresponding to the
transformation:
bar.py 300 bar01.lmp bar10.lmp
The outputs are in kcal/mol and can be compared with the experimental
value of 1.2 kcal/mol and with a simulation value from the literature
(using a different force field): 0.8 kcal/mol
[Gough, Pearlman, Kolmann, J Chem Phys 99 (1993) 9103]. This small
free energy difference is difficult to predict.
These example calculations are for tutorial purposes only. The results
may not be of research quality (not enough sampling, size of the step
in lambda or of the delta for numerical derivative not optimized, no
evaluation of ideal-gas contributions, etc.)

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,89 @@
# created by fftool
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
special_bonds lj/coul 0.0 0.0 0.5
read_data data.0.lmp
# read_restart restart.*.lmp
# reset_timestep 0
pair_style hybrid lj/cut/coul/long 10.0 10.0 lj/cut/coul/long/soft 1 0.5 10.0 10.0 10.0
pair_modify tail yes
kspace_style pppm 1.0e-4
pair_coeff 1 1 none # C-CF C-CF
pair_coeff 1 2 none # C-CF H-D
pair_coeff 1 3 none # C-CF FD-
pair_coeff 1 4 lj/cut/coul/long 0.0000 1.8708 # C-CF Hw
pair_coeff 1 5 lj/cut/coul/long 0.1013 3.3286 # C-CF Ow
pair_coeff 2 2 none # H-D H-D
pair_coeff 2 3 none # H-D FD-
pair_coeff 2 4 lj/cut/coul/long/soft 0.0000 1.5811 1.0 # H-D Hw
pair_coeff 2 5 lj/cut/coul/long/soft 0.0683 2.8131 1.0 # H-D Ow
pair_coeff 3 3 none # FD- FD-
pair_coeff 3 4 lj/cut/coul/long/soft 0.0000 1.7176 0.0 # FD- Hw
pair_coeff 3 5 lj/cut/coul/long/soft 0.0908 3.0559 0.0 # FD- Ow
pair_coeff 4 4 lj/cut/coul/long 0.0000 1.0000 # Hw Hw
pair_coeff 4 5 lj/cut/coul/long 0.0000 1.7792 # Hw Ow
pair_coeff 5 5 lj/cut/coul/long 0.1554 3.1655 # Ow Ow
variable nsteps equal 500000
variable nprint equal ${nsteps}/500
variable ndump equal ${nsteps}/100
# variable nrestart equal ${nsteps}/10
variable temp equal 300.0
variable press equal 1.0
neighbor 2.0 bin
timestep 2.0
velocity all create ${temp} 12345
thermo_style multi
thermo ${nprint}
fix fSHAKE all shake 0.0001 20 ${nprint} b 3 a 4
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 500
run 200000
reset_timestep 0
variable dlambda equal 1.0
variable minusdl equal -1.0
variable dqC equal 0.48*v_dlambda-0.24*(1.0-v_dlambda)
variable dqH equal 0.06*(1.0-v_dlambda)
variable dqF equal -0.12*v_dlambda
compute cFEP all fep ${temp} &
pair lj/cut/coul/long/soft lambda 2 4*5 v_minusdl &
pair lj/cut/coul/long/soft lambda 3 4*5 v_dlambda &
atom charge 1 v_dqC &
atom charge 2 v_dqH &
atom charge 3 v_dqF
fix fFEP all ave/time 1 1 100 c_cFEP[1] c_cFEP[2] file bar01.lmp
# compute cRDF all rdf 100 1 1
# fix fRDF all ave/time 20 100 ${nsteps} c_cRDF file rdf.lammps mode vector
# compute cMSD all msd
# fix fMSD all ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C H F H O
# restart ${nrestart} restart.*.lmp
run ${nsteps}
# write_restart restart.*.lmp
write_data data.*.lmp

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,89 @@
# created by fftool
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
special_bonds lj/coul 0.0 0.0 0.5
read_data data.1.lmp
# read_restart restart.*.lmp
# reset_timestep 0
pair_style hybrid lj/cut/coul/long 10.0 10.0 lj/cut/coul/long/soft 1 0.5 10.0 10.0 10.0
pair_modify tail yes
kspace_style pppm 1.0e-4
pair_coeff 1 1 none # CF-C CF-C
pair_coeff 1 2 none # CF-C HD-
pair_coeff 1 3 none # CF-C F-D
pair_coeff 1 4 lj/cut/coul/long 0.0000 1.8708 # CF-C Hw
pair_coeff 1 5 lj/cut/coul/long 0.1013 3.3286 # CF-C Ow
pair_coeff 2 2 none # HD- HD-
pair_coeff 2 3 none # HD- F-D
pair_coeff 2 4 lj/cut/coul/long/soft 0.0000 1.5811 0.0 # HD- Hw
pair_coeff 2 5 lj/cut/coul/long/soft 0.0683 2.8131 0.0 # HD- Ow
pair_coeff 3 3 none # F-D F-D
pair_coeff 3 4 lj/cut/coul/long/soft 0.0000 1.7176 1.0 # F-D Hw
pair_coeff 3 5 lj/cut/coul/long/soft 0.0908 3.0559 1.0 # F-D Ow
pair_coeff 4 4 lj/cut/coul/long 0.0000 1.0000 # Hw Hw
pair_coeff 4 5 lj/cut/coul/long 0.0000 1.7792 # Hw Ow
pair_coeff 5 5 lj/cut/coul/long 0.1554 3.1655 # Ow Ow
variable nsteps equal 500000
variable nprint equal ${nsteps}/500
variable ndump equal ${nsteps}/100
# variable nrestart equal ${nsteps}/10
variable temp equal 300.0
variable press equal 1.0
neighbor 2.0 bin
timestep 2.0
velocity all create ${temp} 12345
thermo_style multi
thermo ${nprint}
fix fSHAKE all shake 0.0001 20 ${nprint} b 3 a 4
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 500
run 200000
reset_timestep 0
variable dlambda equal 1.0
variable minusdl equal -1.0
variable dqC equal -0.24*v_dlambda+0.48*(1.0-v_dlambda)
variable dqH equal 0.06*v_dlambda
variable dqF equal -0.12*(1.0-v_dlambda)
compute cFEP all fep ${temp} &
pair lj/cut/coul/long/soft lambda 2 4*5 v_dlambda &
pair lj/cut/coul/long/soft lambda 3 4*5 v_minusdl &
atom charge 1 v_dqC &
atom charge 2 v_dqH &
atom charge 3 v_dqF
fix fFEP all ave/time 1 1 100 c_cFEP[1] c_cFEP[2] file bar10.lmp
# compute cRDF all rdf 100 1 1
# fix fRDF all ave/time 20 100 ${nsteps} c_cRDF file rdf.lammps mode vector
# compute cMSD all msd
# fix fMSD all ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C H F H O
# restart ${nrestart} restart.*.lmp
run ${nsteps}
# write_restart restart.*.lmp
write_data data.*.lmp

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,16 @@
CF4->CH4
1 CF-C
2 HD- 1 rCH
3 HD- 1 rCH 2 109.5
4 HD- 1 rCH 2 109.5 3 120.0
5 HD- 1 rCH 2 109.5 3 -120.0
6 F-D 1 rCF 2 109.5 3 -120.0
7 F-D 1 rCF 2 109.5 3 120.0
8 F-D 1 rCF 7 109.5 5 120.0
9 F-D 1 rCF 7 109.5 5 -120.0
rCH = 1.090
rCF = 1.332
CH4-CF4.ff

View File

@ -0,0 +1,28 @@
# oplsaa.ff (kJ/mol, A, deg) version 2013/12/07
ATOMS
# type m/uma q/e pot pars
# alkanes JACS 118(1996)11225, JPC 100(1996)18010
C-CF CT 12.011 -0.24 lj 3.50 0.27614
H-D HC 6.008 0.06 lj 2.50 0.12552
HD- HC 6.008 0.00 lj 2.50 0.12552
# perfluoroalkanes JPCA105:4118(2001)
CF-C CT 12.011 0.48 lj 3.50 0.27614
FD- F 18.998 0.00 lj 2.95 0.22175
F-D F 18.998 -0.12 lj 2.95 0.22175
BONDS
# i j pot re/A ka/kJmol-1
# alkanes OPLS-AA JACS 118(1996)11225, JPC 100(1996)18010
HC CT harm 1.090 2845.12
# perfluoroalkanes JPCA105(2001)4118
F CT harm 1.332 3071.1
ANGLES
# i j k pot th/deg ka/kjmol-1
# alkanes JACS 118(1996)11225, JPC 100(1996)18010
HC CT HC harm 107.80 276.14
# perfluorocarbons JPCA105(2001)4118
F CT F harm 109.1 644.3
# special
HC CT F harm 0.0 0.0

View File

@ -0,0 +1,16 @@
CH4->CF4
1 C-CF
2 H-D 1 rCH
3 H-D 1 rCH 2 109.5
4 H-D 1 rCH 2 109.5 3 120.0
5 H-D 1 rCH 2 109.5 3 -120.0
6 FD- 1 rCF 2 109.5 3 -120.0
7 FD- 1 rCF 2 109.5 3 120.0
8 FD- 1 rCF 7 109.5 5 120.0
9 FD- 1 rCF 7 109.5 5 -120.0
rCH = 1.090
rCF = 1.332
CH4-CF4.ff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,14 @@
# SPC/E water (kJ/mol, A, deg)
ATOMS
# type m/uma q/e pot pars
Ow Ow 15.999 -0.8476 lj 3.1655 0.6503
Hw Hw 1.008 0.4238 lj 1.0000 0.0000
BONDS
# i j pot re/A ka/kJmol-1
Ow Hw cons 1.000 4331.53
ANGLES
# i j k pot th/deg ka/kjmol-1
Hw Ow Hw cons 109.47 317.57

View File

@ -0,0 +1,7 @@
SPCE
Hw
Ow 1 1.000
Hw 2 1.000 1 109.47
spce.ff

View File

@ -0,0 +1,59 @@
Free Energy of Hydration of Methane
===================================
Example calculation of the free energy of hydration of methane with
LAMMPS using *compute fep* and *fix adapt/fep*.
Methane is represented by the 5-site OPLS-AA model (1 molecule). Water
is represented by the 4-site TIP4P-Ew model (360 molecules). Interactions
of sites that are being created or deleted are treated using soft-core
verions of the Lennard-Jones and Coulomb potentials (*pair
lj/cut/coul/long/soft*) in order to avoid singularities.
The following directories contain input files and results for
calculations using free-energy perturbation (FEP) and
finite-difference thermodynamic integration (FDTI):
* `mols` -- molecule description files and force field database used
to create the initial configuration used for the simulations
`data.lmp`
* `fep01` -- Calculation using FEP, multi-stage creation of a methane
molecule. Results in `fep01.lmp`
* `fep10` -- Calculation using FEP, multi-stage deletion of a methane
molecule. Results in `fep10.lmp`
* `fdti01` -- Calculation using TI/FDTI, creation of a methane
molecule. Results in `fdti01.lmp`
* `fdti10` -- Calculation using TI/FDTI, deletion a methane
molecule. Results in `fdti10.lmp`
The free-energy profiles can be observed by plotting the values in the
third column of the results files. The Python scripts `fep.py`,
`nti.py` and `fdti.py` found in the `tools` directory can be used to
calculate the free-energy differences corresponding to the above
transformations:
fep.py 300 < fep01.lmp
fep.py 300 < fep10.lmp
nti.py 300 0.002 < fdti01.lmp
nti.py 300 0.002 < fdti10.lmp
fdti.py 300 0.002 < fdti01.lmp
fdti.py 300 0.002 < fdti10.lmp
The outputs are in kcal/mol and can be compared with the experimental
value of 2.0 kcal/mol, or with a simulation value from the literature
obtained with the same force field models used here: 2.27 kcal/mol
[MR Shirts, VS Pande, J Chem Phys 122 (2005) 134508].
These example calculations are for tutorial purposes only. The results
may not be of research quality (not enough sampling, size of the step
in lambda or of the delta for numerical derivative not optimized, no
evaluation of ideal-gas contributions, etc.)

View File

@ -0,0 +1,22 @@
# Time-averaged data for fix fFEP
# TimeStep c_cFEP[1] c_cFEP[2]
125000 0.000490778 0.999186
250000 0.00619786 0.989679
375000 0.0128045 0.978813
500000 0.0203102 0.96667
625000 0.0272161 0.955808
750000 0.0215735 0.965087
875000 0.0148509 0.975821
1000000 0.00936378 0.984752
1125000 0.00490088 0.992007
1250000 0.00388501 0.99366
1375000 0.00122887 0.998043
1500000 0.000498063 0.999247
1625000 -0.000468616 1.00085
1750000 -0.00131731 1.00226
1875000 -0.00247035 1.00419
2000000 -0.00306799 1.00519
2125000 -0.00355565 1.006
2250000 -0.00434289 1.00733
2375000 -0.00497634 1.00839
2500000 -0.00576373 1.00972

View File

@ -0,0 +1,80 @@
# created by fftool
print "Methane in TIP4P water"
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
special_bonds lj/coul 0.0 0.0 0.5
read_data data.lmp
pair_style hybrid &
lj/cut/coul/long 10.0 10.0 &
lj/cut/tip4p/long/soft 3 4 2 2 0.125 1 0.5 10.0 10.0 10.0
pair_modify tail no
kspace_style pppm/tip4p 1.0e-4
pair_coeff 1 1 lj/cut/coul/long 0.0660 3.5000 # C4H C4H
pair_coeff 1 2 lj/cut/coul/long 0.0445 2.9580 # C4H H
pair_coeff 2 2 lj/cut/coul/long 0.0300 2.5000 # H H
pair_coeff 3 3 lj/cut/tip4p/long/soft 0.1628 3.1644 1.0 # Ow Ow
pair_coeff 3 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # Ow Hw
pair_coeff 4 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # Hw Hw
pair_coeff 1 3 lj/cut/tip4p/long/soft 0.1036 3.3279 0.0 # C4H Ow
pair_coeff 1 4 lj/cut/tip4p/long/soft 0.0000 1.0000 0.0 # C4H Hw
pair_coeff 2 3 lj/cut/tip4p/long/soft 0.0699 2.8126 0.0 # H Ow
pair_coeff 2 4 lj/cut/tip4p/long/soft 0.0000 1.0000 0.0 # H Hw
variable nsteps equal 2500000
variable nprint equal ${nsteps}/500
variable ndump equal ${nsteps}/100
variable temp equal 300.0
variable press equal 1.0
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 a 2
neighbor 2.0 bin
timestep 2.0
velocity all create ${temp} 12345
thermo_style multi
thermo ${nprint}
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 1000
run 250000
reset_timestep 0
variable lambda equal ramp(0.0,1.0)
fix fADAPT all adapt/fep 125000 pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_lambda after yes
variable dlambda equal 0.002
compute cFEP all fep ${temp} pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_dlambda
fix fFEP all ave/time 25 4000 125000 c_cFEP[1] c_cFEP[2] file fdti01.lmp
# compute cRDF all rdf 100 3 3 3 4 4 4
# fix fRDF all ave/time 2000 500 ${nsteps} c_cRDF file rdf.lammps mode vector
group owater type 3
compute cMSD owater msd
fix fMSD owater ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C H O H
run ${nsteps}
# write_restart restart.*.lmp
write_data data.*.lmp

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,22 @@
# Time-averaged data for fix fFEP
# TimeStep c_cFEP[1] c_cFEP[2]
125000 0.00652155 0.989123
250000 0.00572718 0.990445
375000 0.00506788 0.991544
500000 0.00448303 0.992522
625000 0.00370796 0.993818
750000 0.00277535 0.995385
875000 0.00230833 0.996171
1000000 0.000982551 0.998409
1125000 0.000314678 0.99954
1250000 -0.000839025 1.0015
1375000 -0.00155615 1.00273
1500000 -0.00320598 1.00553
1625000 -0.00587133 1.01011
1750000 -0.00867306 1.01495
1875000 -0.0134768 1.02327
2000000 -0.0229314 1.03995
2125000 -0.0282333 1.04896
2250000 -0.0202762 1.03477
2375000 -0.0123279 1.02096
2500000 -0.00546203 1.00923

View File

@ -0,0 +1,80 @@
# created by fftool
print "Methane in TIP4P water"
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
special_bonds lj/coul 0.0 0.0 0.5
read_data data.lmp
pair_style hybrid &
lj/cut/coul/long 10.0 10.0 &
lj/cut/tip4p/long/soft 3 4 2 2 0.125 1 0.5 10.0 10.0 10.0
pair_modify tail no
kspace_style pppm/tip4p 1.0e-4
pair_coeff 1 1 lj/cut/coul/long 0.0660 3.5000 # C4H C4H
pair_coeff 1 2 lj/cut/coul/long 0.0445 2.9580 # C4H H
pair_coeff 2 2 lj/cut/coul/long 0.0300 2.5000 # H H
pair_coeff 3 3 lj/cut/tip4p/long/soft 0.1628 3.1644 1.0 # Ow Ow
pair_coeff 3 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # Ow Hw
pair_coeff 4 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # Hw Hw
pair_coeff 1 3 lj/cut/tip4p/long/soft 0.1036 3.3279 1.0 # C4H Ow
pair_coeff 1 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # C4H Hw
pair_coeff 2 3 lj/cut/tip4p/long/soft 0.0699 2.8126 1.0 # H Ow
pair_coeff 2 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # H Hw
variable nsteps equal 2500000
variable nprint equal ${nsteps}/500
variable ndump equal ${nsteps}/100
variable temp equal 300.0
variable press equal 1.0
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 a 2
neighbor 2.0 bin
timestep 2.0
velocity all create ${temp} 12345
thermo_style multi
thermo ${nprint}
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 1000
run 250000
reset_timestep 0
variable lambda equal ramp(1.0,0.0)
fix fADAPT all adapt/fep 125000 pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_lambda after yes
variable dlambda equal -0.002
compute cFEP all fep ${temp} pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_dlambda
fix fFEP all ave/time 25 4000 125000 c_cFEP[1] c_cFEP[2] file fdti10.lmp
# compute cRDF all rdf 100 3 3 3 4 4 4
# fix fRDF all ave/time 2000 500 ${nsteps} c_cRDF file rdf.lammps mode vector
group owater type 3
compute cMSD owater msd
fix fMSD owater ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C H O H
run ${nsteps}
# write_restart restart.*.lmp
write_data data.*.lmp

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,22 @@
# Time-averaged data for fix fFEP
# TimeStep c_cFEP[1] c_cFEP[2]
125000 0.0806308 0.881044
250000 0.251798 0.66945
375000 0.448135 0.502774
500000 0.670275 0.384034
625000 0.873864 0.343148
750000 0.689533 0.484145
875000 0.474153 0.600305
1000000 0.303566 0.740249
1125000 0.165693 0.855143
1250000 0.130949 0.88524
1375000 0.0518948 0.976114
1500000 0.0283335 1.0025
1625000 -0.000790228 1.04259
1750000 -0.0263791 1.0736
1875000 -0.0587947 1.12568
2000000 -0.0764999 1.15472
2125000 -0.091264 1.17869
2250000 -0.11294 1.21778
2375000 -0.130721 1.25162
2500000 -0.151808 1.29392

View File

@ -0,0 +1,80 @@
# created by fftool
print "Methane in TIP4P water"
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
special_bonds lj/coul 0.0 0.0 0.5
read_data data.lmp
pair_style hybrid &
lj/cut/coul/long 10.0 10.0 &
lj/cut/tip4p/long/soft 3 4 2 2 0.125 1 0.5 10.0 10.0 10.0
pair_modify tail no
kspace_style pppm/tip4p 1.0e-4
pair_coeff 1 1 lj/cut/coul/long 0.0660 3.5000 # C4H C4H
pair_coeff 1 2 lj/cut/coul/long 0.0445 2.9580 # C4H H
pair_coeff 2 2 lj/cut/coul/long 0.0300 2.5000 # H H
pair_coeff 3 3 lj/cut/tip4p/long/soft 0.1628 3.1644 1.0 # Ow Ow
pair_coeff 3 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # Ow Hw
pair_coeff 4 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # Hw Hw
pair_coeff 1 3 lj/cut/tip4p/long/soft 0.1036 3.3279 0.0 # C4H Ow
pair_coeff 1 4 lj/cut/tip4p/long/soft 0.0000 1.0000 0.0 # C4H Hw
pair_coeff 2 3 lj/cut/tip4p/long/soft 0.0699 2.8126 0.0 # H Ow
pair_coeff 2 4 lj/cut/tip4p/long/soft 0.0000 1.0000 0.0 # H Hw
variable nsteps equal 2500000
variable nprint equal ${nsteps}/500
variable ndump equal ${nsteps}/100
variable temp equal 300.0
variable press equal 1.0
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 a 2
neighbor 2.0 bin
timestep 2.0
velocity all create ${temp} 12345
thermo_style multi
thermo ${nprint}
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 1000
run 250000
reset_timestep 0
variable lambda equal ramp(0.0,1.0)
fix fADAPT all adapt/fep 125000 pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_lambda after yes
variable dlambda equal 0.05
compute cFEP all fep ${temp} pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_dlambda
fix fFEP all ave/time 25 4000 125000 c_cFEP[1] c_cFEP[2] file fep01.lmp
# compute cRDF all rdf 100 3 3 3 4 4 4
# fix fRDF all ave/time 2000 500 ${nsteps} c_cRDF file rdf.lammps mode vector
group owater type 3
compute cMSD owater msd
fix fMSD owater ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C H O H
run ${nsteps}
# write_restart restart.*.lmp
write_data data.*.lmp

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,22 @@
# Time-averaged data for fix fFEP
# TimeStep c_cFEP[1] c_cFEP[2]
125000 0.156203 0.771285
250000 0.136448 0.799233
375000 0.12198 0.820426
500000 0.109453 0.840255
625000 0.0921141 0.868102
750000 0.0713015 0.905735
875000 0.0622431 0.922352
1000000 0.0330546 0.981942
1125000 0.0197866 1.01134
1250000 -0.00412826 1.06765
1375000 -0.0174302 1.11343
1500000 -0.0511646 1.18453
1625000 -0.105676 1.36325
1750000 -0.161616 1.56478
1875000 -0.258263 1.93825
2000000 -0.447487 3.0595
2125000 -0.549377 3.01163
2250000 -0.379009 2.01933
2375000 -0.20934 1.45678
2500000 -0.0644056 1.1252

View File

@ -0,0 +1,80 @@
# created by fftool
print "Methane in TIP4P water"
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
special_bonds lj/coul 0.0 0.0 0.5
read_data data.lmp
pair_style hybrid &
lj/cut/coul/long 10.0 10.0 &
lj/cut/tip4p/long/soft 3 4 2 2 0.125 1 0.5 10.0 10.0 10.0
pair_modify tail no
kspace_style pppm/tip4p 1.0e-4
pair_coeff 1 1 lj/cut/coul/long 0.0660 3.5000 # C4H C4H
pair_coeff 1 2 lj/cut/coul/long 0.0445 2.9580 # C4H H
pair_coeff 2 2 lj/cut/coul/long 0.0300 2.5000 # H H
pair_coeff 3 3 lj/cut/tip4p/long/soft 0.1628 3.1644 1.0 # Ow Ow
pair_coeff 3 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # Ow Hw
pair_coeff 4 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # Hw Hw
pair_coeff 1 3 lj/cut/tip4p/long/soft 0.1036 3.3279 1.0 # C4H Ow
pair_coeff 1 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # C4H Hw
pair_coeff 2 3 lj/cut/tip4p/long/soft 0.0699 2.8126 1.0 # H Ow
pair_coeff 2 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # H Hw
variable nsteps equal 2500000
variable nprint equal ${nsteps}/500
variable ndump equal ${nsteps}/100
variable temp equal 300.0
variable press equal 1.0
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 a 2
neighbor 2.0 bin
timestep 2.0
velocity all create ${temp} 12345
thermo_style multi
thermo ${nprint}
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 1000
run 250000
reset_timestep 0
variable lambda equal ramp(1.0,0.0)
fix fADAPT all adapt/fep 125000 pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_lambda after yes
variable dlambda equal -0.05
compute cFEP all fep ${temp} pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_dlambda
fix fFEP all ave/time 25 4000 125000 c_cFEP[1] c_cFEP[2] file fep10.lmp
# compute cRDF all rdf 100 3 3 3 4 4 4
# fix fRDF all ave/time 2000 500 ${nsteps} c_cRDF file rdf.lammps mode vector
group owater type 3
compute cMSD owater msd
fix fMSD owater ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C H O H
run ${nsteps}
# write_restart restart.*.lmp
write_data data.*.lmp

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,9 @@
methane
1 C4H
2 H 1 1.09
3 H 1 1.09 2 109.5
4 H 1 1.09 2 109.5 3 120.0
5 H 1 1.09 2 109.5 3 -120.0
oplsaa.ff

View File

@ -0,0 +1,85 @@
# oplsaa.ff (kJ/mol, A, deg) version 2013/07/27
ATOMS
# type m/uma q/e pot pars
# alkanes OPLS-AA JACS 118(1996)11225, JPC 100(1996)18010
C4C CT 12.011 0.00 lj 3.50 0.27614
C4H CT 4.011 -0.24 lj 3.50 0.27614
C3H CT 12.011 -0.18 lj 3.50 0.27614
C2H CT 12.011 -0.12 lj 3.50 0.27614
H HC 4.008 0.06 lj 2.50 0.12552
# alcohols OPLS-AA JACS 118(1996)11225, JPC 100(1996)18010
CTO CT 12.011 0.145 lj 3.50 0.27614
CET CT 12.011 -0.140 lj 3.50 0.27614
C2T CT 12.011 -0.080 lj 3.50 0.27614
H1O HC 1.008 0.040 lj 2.50 0.12552
OH OH 15.999 -0.683 lj 3.12 0.71128
HO HO 1.008 0.418 lj 0.00 0.00000
# aromatics OPLS-AA JACS 118(1996)11225
CA CA 12.011 -0.115 lj 3.55 0.43932
HA HA 1.008 0.115 lj 2.42 0.12552
# toluene OPLS-AA JACS 118(1996)11217, JCC 25(2004)1322
CAP CA 12.011 -0.154 lj 3.55 0.43932
CAM CA 12.011 -0.145 lj 3.55 0.43932
CAO CA 12.011 -0.147 lj 3.55 0.43932
CAT CA 12.011 -0.081 lj 3.55 0.43932
CTT CT 12.011 -0.197 lj 3.50 0.27614
HAT HA 1.008 0.148 lj 2.42 0.12552
HT HC 1.008 0.092 lj 2.42 0.12552
BONDS
# i j pot re/A ka/kJmol-1
# alkanes OPLS-AA JACS 118(1996)11225, JPC 100(1996)18010
HC CT harm 1.090 2845.12
CT CT harm 1.529 2242.62
# alcohols OPLS-AA JACS 118(1996)11225, JPC 100(1996)18010
CT OH harm 1.410 2677.76
HO OH cons 0.945 4627.50
# aromatics AMBER JACS 117(1995)5179
CA CA harm 1.400 3924.6
CA HA cons 1.080 3071.1
# toluene AMBER JACS 117(1995)5179
CA CT harm 1.510 2652.7
ANGLES
# i j k pot th/deg ka/kjmol-1
# alkanes OPLS-AA JACS 118(1996)11225, JPC 100(1996)18010
HC CT HC harm 107.80 276.14
HC CT CT harm 110.70 313.80
CT CT CT harm 112.70 488.27
# alcohols OPLS-AA JACS 118(1996)11225, JPC 100(1996)18010
CT CT OH harm 109.50 418.40
HC CT OH harm 109.50 292.88
CT OH HO harm 108.50 460.24
# aromatics AMBER JACS 117(1995)5179
CA CA HA harm 120.00 292.88
CA CA CA harm 120.00 527.18
# toluene AMBER JACS 117(1995)5179
CA CA CT harm 120.00 585.76
CA CT HC harm 109.50 418.4
DIHEDRALS
# i j k l pot v1 v2 v3 v4
# alkanes OPLS-AA JACS 118: 11225 (1996); J Phys Chem 100: 18010 (1996)
HC CT CT HC opls 0.0000 0.0000 1.3305 0.0000
HC CT CT CT opls 0.0000 0.0000 1.5313 0.0000
CT CT CT CT opls 7.2800 -0.6569 1.1673 0.0000
# alcohols OPLS-AA JACS 118(1996)11225, AMBER98 (OCCO) 117(1995)5179
HC CT OH HO opls 0.0000 0.0000 1.8828 0.0000
CT CT OH HO opls -1.4895 -0.7280 2.0585 0.0000
HC CT CT OH opls 0.0000 0.0000 1.9581 0.0000
OH CT CT OH opls 0.0000 -9.8324 1.2050 0.0000
# aromatics AMBER JACS 117(1995)5179
CA CA CA CA opls 0.0000 30.334 0.0000 0.0000
CA CA CA HA opls 0.0000 30.334 0.0000 0.0000
HA CA CA HA opls 0.0000 30.334 0.0000 0.0000
# toluene AMBER JACS 117: 5179 (1995)
CT CA CA CA opls 0.0000 30.334 0.0000 0.0000
CT CA CA HA opls 0.0000 30.334 0.0000 0.0000
CA CA CT HC opls 0.0000 0.0000 0.0000 0.0000
IMPROPER
# improper aromatics AMBER JACS 117(1995)5179
CA CA CA HA opls 0.0000 9.2048 0.0000 0.0000
# improper toluene AMBER JACS 117: 5179 (1995)
CA CA CA CT opls 0.0000 9.2048 0.0000 0.0000

View File

@ -0,0 +1,14 @@
# TIP4P-Ew water (kJ/mol, A, deg)
ATOMS
# type m/uma q/e pot pars
Ow Ow 15.999 -1.04844 lj 3.16435 0.680946
Hw Hw 1.008 0.52422 lj 0.00000 0.000000
BONDS
# i j pot re/A ka/kJmol-1
Ow Hw cons 0.9572 4331.53
ANGLES
# i j k pot th/deg ka/kjmol-1
Hw Ow Hw cons 104.52 317.57

View File

@ -0,0 +1,7 @@
TIP4P-Ew
Ow
Hw 1 0.9572
Hw 1 0.9572 2 104.52
tip4pew.ff

11
examples/USER/fep/README Normal file
View File

@ -0,0 +1,11 @@
compute fep examples
====================
* `CH4hyd` -- free energy of hydration of methane (very simple). FEP
and TI methods.
* `CH4-CF4` -- free energy difference of transforming methane into
perfluoromethane, in water (simple). BAR method.
* `CC-CO` -- free energy difference of transforming ethane into
methanol, in water (more complex). FEP, TI and BAR methods.