Added first go at checks on PLUMED interface to LAMMPS

examples/USER/plumed/check-plumed.sh

@@ -1,5 +1,10 @@
 #!/bin/bash
 
+LAMMPS=../../../src/lmp_mpi
+
+# Run first LAMMPS calculation
+$LAMMPS < in.peptide-plumed
+
 # Check PLUMED positions
 nlines=`paste plumed.xyz lammps.xyz | awk '{if( $2<$6-0.0001 || $2>$6+0.0001 || $3<$7-0.0001 || $3>$7+0.0001 || $4<$8-0.0002 || $4>$8+0.0002 ) if( $5!="Timestep:" && $1!=2004 ) print $0}' | wc -l` 
 if [ "$nlines" -gt 0 ] ; then
@@ -30,3 +35,47 @@ if [ "$nlines" -gt 0 ] ; then
 fi
 rm -f plmd_energy2
 
+# Check PLMD mass and charge
+nlines=`wc -l mq_lammps | awk '{print $1}'`
+sline=`grep -n "mass q" mq_lammps | awk '{print $1}' | sed -e s/:ITEM://`
+for ((i=$sline+1;i<$nlines;i++)); do
+    # Mass and charge from LAMMPS
+    index=`head -n $i mq_lammps | tail -n 1 | awk '{print $1}'`
+    l_mass=`head -n $i mq_lammps | tail -n 1 | awk '{print $2}'`
+    l_charge=`head -n $i mq_lammps | tail -n 1 | awk '{print $3}'`
+    # Mass and charge from PLUMED
+    p_mass=`head -n $(($index+1)) mq_plumed | tail -n 1 | awk '{print $2}'`
+    p_charge=`head -n $(($index+1)) mq_plumed | tail -n 1 | awk '{print $3}'`
+    # Check PLUMED mass is same as lammps mass
+    mdiff=`echo \( $l_mass - $p_mass \) \> 0 | bc -l` 
+    if [ "$mdiff" -gt 0 ] ; then
+         echo ERROR passing masses from LAMMPS to PLUMED
+    fi 
+    # Check PLUMED charge is same as lammps charge
+    qdiff=`echo \( $l_charge - $p_charge \) \> 0 | bc -l` 
+    if [ "$qdiff" -gt 0 ] ; then
+         echo ERROR passing charges from LAMMPS to PLUMED
+    fi 
+done
+
+# Run calculations to test adding restraint on bond
+$LAMMPS < in.peptide-plumed-plumed-restraint
+$LAMMPS < in.peptide-plumed-lammps-restraint
+# Now compare value of distance when lammps and plumed restraint the distance
+nlines=`paste lammps_restraint plumed_restraint | tail -n +2 | awk '{if( $2<$4-0.0001 || $2>$4+0.0001 ) print $0}' | wc -l`
+if [ "$nlines" -gt 0 ] ; then
+      echo ERROR passing forces from PLUMED back to LAMMPS
+fi
+
+# Nothing from here works
+
+# Now try to simply increase the size of the box by applying a moving restraint on the volume
+$LAMMPS < in.peptide-plumed-expand
+
+# Now run calculations to test virial
+$LAMMPS < in.peptide-plumed-npt
+$LAMMPS < in.peptide-plumed-npt2	
+
+# Now run calculations to check forces on energy
+$LAMMPS < in.peptide-plumed-engforce-ref
+$LAMMPS < in.peptide-plumed-eng-force-plumed

examples/USER/plumed/clear-files.sh

@@ -0,0 +1,8 @@
+#!/bin/bash
+
+# Data from first set of checks
+rm bck.* plmd_energy lammps_energy mq_plumed mq_lammps lammps.xyz plumed.xyz p.log
+# Data from checks on restraints
+rm bck.* p.log lammps_restraint plumed_restraint
+# Data from checks on virial
+rm bck.* lammps_energy lammps.xyz log.lammps plmd_volume p.log plmd_volume_without_restraint plmd_volume_with_restraint

examples/USER/plumed/in.peptide-plumed

@@ -39,5 +39,6 @@ dump            dd all xyz 10 lammps.xyz
 variable        step equal step
 variable        pe equal pe
 fix             5 all print 10 "$(v_step) $(v_pe)" file lammps_energy
+dump            mq all custom 200 mq_lammps id mass q 
 
 run		101

examples/USER/plumed/in.peptide-plumed-eng-force-plumed

@@ -0,0 +1,41 @@
+# Solvated 5-mer peptide
+
+units		real
+atom_style	full
+
+pair_style	lj/charmm/coul/long 8.0 10.0 10.0
+bond_style      harmonic
+angle_style     charmm
+dihedral_style  charmm
+improper_style  harmonic
+kspace_style	pppm 0.0001
+
+read_data	data.peptide
+
+neighbor	2.0 bin
+neigh_modify	delay 5
+
+timestep	1.8181818181818181
+
+group		peptide type <= 12
+group		one id 2 4 5 6
+group		two id 80 82 83 84
+group		ref id 37
+group		colvar union one two ref
+
+fix		1 all nvt temp  363.0 363.0 90.90909090909091 tchain 1
+
+fix		2 all plumed plumedfile plumed-engforce.dat outfile p.log
+
+#dump		1 colvar custom 1 dump.colvar.lammpstrj id xu yu zu fx fy fz
+#dump_modify 1 sort id
+
+thermo_style	custom step temp etotal pe ke epair ebond f_2
+thermo		10
+dump            dd all xyz 10 lammps.xyz
+variable        step equal step
+variable        pe equal pe
+fix             5 all print 10 "$(v_step) $(v_pe)" file lammps_energy
+dump            mq all custom 200 mq_lammps id mass q 
+
+run		101

examples/USER/plumed/in.peptide-plumed-engforce-ref

@@ -0,0 +1,40 @@
+# Solvated 5-mer peptide
+
+units		real
+atom_style	full
+
+pair_style	lj/charmm/coul/long 8.0 10.0 10.0
+bond_style      harmonic
+angle_style     charmm
+dihedral_style  charmm
+improper_style  harmonic
+kspace_style	pppm 0.0001
+
+read_data	data.peptide
+
+neighbor	2.0 bin
+neigh_modify	delay 5
+
+timestep	2.0
+
+group		peptide type <= 12
+group		one id 2 4 5 6
+group		two id 80 82 83 84
+group		ref id 37
+group		colvar union one two ref
+
+fix		1 all nvt temp  300.0 300.0 100.0 tchain 1
+
+fix		2 all plumed plumedfile plumed-eng-ref.dat outfile p.log
+
+#dump		1 colvar custom 1 dump.colvar.lammpstrj id xu yu zu fx fy fz
+#dump_modify 1 sort id
+
+thermo_style	custom step temp etotal pe ke epair ebond f_2
+thermo		10
+dump            dd all xyz 10 lammps.xyz
+variable        step equal step
+variable        pe equal pe
+dump            mq all custom 200 mq_lammps id mass q 
+
+run		101

examples/USER/plumed/in.peptide-plumed-expand

@@ -0,0 +1,40 @@
+# Solvated 5-mer peptide
+
+units		real
+atom_style	full
+
+pair_style	lj/charmm/coul/long 8.0 10.0 10.0
+bond_style      harmonic
+angle_style     charmm
+dihedral_style  charmm
+improper_style  harmonic
+kspace_style	pppm 0.0001
+
+read_data	data.peptide
+
+neighbor	2.0 bin
+neigh_modify	delay 5
+
+timestep	2.0
+
+group		peptide type <= 12
+group		one id 2 4 5 6
+group		two id 80 82 83 84
+group		ref id 37
+group		colvar union one two ref
+
+fix		1 all npt temp  275.0 275.0 100.0 iso 0.987 0.987 400.0 tchain 1
+
+fix		2 all plumed plumedfile plumed-expand.dat outfile p.log
+
+#dump		1 colvar custom 1 dump.colvar.lammpstrj id xu yu zu fx fy fz
+#dump_modify 1 sort id
+
+thermo_style	custom step temp etotal pe ke epair ebond f_2
+thermo		10
+dump            dd all xyz 10 lammps.xyz
+variable        step equal step
+variable        pe equal pe
+dump            mq all custom 200 mq_lammps id mass q 
+
+run		101

examples/USER/plumed/in.peptide-plumed-lammps-restraint

@@ -0,0 +1,44 @@
+# Solvated 5-mer peptide
+
+units		real
+atom_style	full
+
+pair_style	lj/charmm/coul/long 8.0 10.0 10.0
+bond_style      harmonic
+angle_style     charmm
+dihedral_style  charmm
+improper_style  harmonic
+kspace_style	pppm 0.0001
+
+read_data	data.peptide
+
+neighbor	2.0 bin
+neigh_modify	delay 5
+
+timestep	2.0
+
+group		peptide type <= 12
+group		one id 2 4 5 6
+group		two id 80 82 83 84
+group		ref id 37
+group		colvar union one two ref
+
+fix		1 all nvt temp  275.0 275.0 100.0 tchain 1
+fix             6 all restrain bond 45 48 1000.0 1000.0 6.0
+fix		2 all plumed plumedfile plumed-norestraint.dat outfile p.log
+fix		2a ref setforce 0.0 0.0 0.0
+
+fix		4 all shake 0.0001 10 100 b 4 6 8 10 12 14 18 a 31
+
+#dump		1 colvar custom 1 dump.colvar.lammpstrj id xu yu zu fx fy fz
+#dump_modify 1 sort id
+
+thermo_style	custom step temp etotal pe ke epair ebond f_2
+thermo		10
+dump            dd all xyz 10 lammps.xyz
+variable        step equal step
+variable        pe equal pe
+fix             5 all print 10 "$(v_step) $(v_pe)" file lammps_energy
+dump            mq all custom 200 mq_lammps id mass q 
+
+run		101

examples/USER/plumed/in.peptide-plumed-npt

@@ -0,0 +1,43 @@
+# Solvated 5-mer peptide
+
+units		real
+atom_style	full
+
+pair_style	lj/charmm/coul/long 8.0 10.0 10.0
+bond_style      harmonic
+angle_style     charmm
+dihedral_style  charmm
+improper_style  harmonic
+kspace_style	pppm 0.0001
+
+read_data	data.peptide
+
+neighbor	2.0 bin
+neigh_modify	delay 5
+
+timestep	2.0
+
+group		peptide type <= 12
+group		one id 2 4 5 6
+group		two id 80 82 83 84
+group		ref id 37
+group		colvar union one two ref
+
+fix             4 all shake 0.0001 10 100 b 4 6 8 10 12 14 18 a 31
+fix		1 all npt temp  275.0 275.0 100.0 iso 0.987 0.987 400.0 tchain 1
+
+fix		2 all plumed plumedfile plumed_npt.dat outfile p.log
+fix		2a ref setforce 0.0 0.0 0.0
+
+#dump		1 colvar custom 1 dump.colvar.lammpstrj id xu yu zu fx fy fz
+#dump_modify 1 sort id
+
+thermo_style	custom step temp etotal pe ke epair ebond f_2
+thermo		10
+dump            dd all xyz 10 lammps.xyz
+variable        step equal step
+variable        pe equal pe
+fix             5 all print 10 "$(v_step) $(v_pe)" file lammps_energy
+dump            mq all custom 200 mq_lammps id mass q 
+
+run		101

examples/USER/plumed/in.peptide-plumed-npt2

@@ -0,0 +1,43 @@
+# Solvated 5-mer peptide
+
+units		real
+atom_style	full
+
+pair_style	lj/charmm/coul/long 8.0 10.0 10.0
+bond_style      harmonic
+angle_style     charmm
+dihedral_style  charmm
+improper_style  harmonic
+kspace_style	pppm 0.0001
+
+read_data	data.peptide
+
+neighbor	2.0 bin
+neigh_modify	delay 5
+
+timestep	2.0
+
+group		peptide type <= 12
+group		one id 2 4 5 6
+group		two id 80 82 83 84
+group		ref id 37
+group		colvar union one two ref
+
+fix             4 all shake 0.0001 10 100 b 4 6 8 10 12 14 18 a 31
+fix		1 all npt temp  275.0 275.0 100.0 iso 987.9 987.9 400.0 tchain 1
+
+fix		2 all plumed plumedfile plumed_npt2.dat outfile p.log
+fix		2a ref setforce 0.0 0.0 0.0
+
+#dump		1 colvar custom 1 dump.colvar.lammpstrj id xu yu zu fx fy fz
+#dump_modify 1 sort id
+
+thermo_style	custom step temp etotal pe ke epair ebond f_2
+thermo		10
+dump            dd all xyz 10 lammps.xyz
+variable        step equal step
+variable        pe equal pe
+fix             5 all print 10 "$(v_step) $(v_pe)" file lammps_energy
+dump            mq all custom 200 mq_lammps id mass q 
+
+run		101

examples/USER/plumed/in.peptide-plumed-plumed-restraint

@@ -0,0 +1,44 @@
+# Solvated 5-mer peptide
+
+units		real
+atom_style	full
+
+pair_style	lj/charmm/coul/long 8.0 10.0 10.0
+bond_style      harmonic
+angle_style     charmm
+dihedral_style  charmm
+improper_style  harmonic
+kspace_style	pppm 0.0001
+
+read_data	data.peptide
+
+neighbor	2.0 bin
+neigh_modify	delay 5
+
+timestep	2.0
+
+group		peptide type <= 12
+group		one id 2 4 5 6
+group		two id 80 82 83 84
+group		ref id 37
+group		colvar union one two ref
+
+fix		1 all nvt temp  275.0 275.0 100.0 tchain 1
+
+fix		2 all plumed plumedfile plumed-restraint.dat outfile p.log
+fix		2a ref setforce 0.0 0.0 0.0
+
+fix		4 all shake 0.0001 10 100 b 4 6 8 10 12 14 18 a 31
+
+#dump		1 colvar custom 1 dump.colvar.lammpstrj id xu yu zu fx fy fz
+#dump_modify 1 sort id
+
+thermo_style	custom step temp etotal pe ke epair ebond f_2
+thermo		10
+dump            dd all xyz 10 lammps.xyz
+variable        step equal step
+variable        pe equal pe
+fix             5 all print 10 "$(v_step) $(v_pe)" file lammps_energy
+dump            mq all custom 200 mq_lammps id mass q 
+
+run		101

examples/USER/plumed/plumed-eng-ref.dat

@@ -0,0 +1,2 @@
+e: ENERGY
+PRINT ARG=e FILE=energy_ref FMT=%8.4f

examples/USER/plumed/plumed-engforce.dat

@@ -0,0 +1,3 @@
+e: ENERGY
+RESTRAINT ARG=e AT=0 SLOPE=0.1
+PRINT ARG=e FILE=engforced FMT=%8.4f

examples/USER/plumed/plumed-expand.dat

@@ -0,0 +1,3 @@
+v: VOLUME 
+r: MOVINGRESTRAINT ARG=v KAPPA0=100. AT0=20 AT1=30 STEP0=0 STEP1=100
+PRINT ARG=v,r.* FILE=expanding

examples/USER/plumed/plumed-norestraint.dat

@@ -0,0 +1,3 @@
+UNITS LENGTH=A
+dd: DISTANCE ATOMS=45,48 
+PRINT ARG=dd FILE=lammps_restraint FMT=%8.4f

examples/USER/plumed/plumed-restraint.dat

@@ -0,0 +1,4 @@
+UNITS LENGTH=A ENERGY=kcal/mol
+dd: DISTANCE ATOMS=45,48 
+RESTRAINT ARG=dd KAPPA=2000 AT=6.0
+PRINT ARG=dd FILE=plumed_restraint FMT=%8.4f 

examples/USER/plumed/plumed.dat

@@ -4,3 +4,4 @@ c2: COM ATOMS=80,82,83,84
 energy: ENERGY
 DUMPATOMS ATOMS=1-2004 FILE=plumed.xyz UNITS=A PRECISION=4 STRIDE=10
 PRINT ARG=energy STRIDE=10 FMT=%8.4f FILE=plmd_energy
+DUMPMASSCHARGE FILE=mq_plumed 

examples/USER/plumed/plumed_npt.dat

@@ -0,0 +1,2 @@
+vv: VOLUME
+PRINT ARG=vv FILE=plmd_volume_without_restraint

examples/USER/plumed/plumed_npt2.dat

@@ -0,0 +1,3 @@
+vv: VOLUME
+RESTRAINT AT=0.0 ARG=vv SLOPE=-60.221429
+PRINT ARG=vv FILE=plmd_volume_with_restraint