Tidied up example directory for PLUMED

#  Needs make yes-molecule

#        make yes-rigid

#        make yes-kspace

#        make yes-user-plumed    ( obviously )

#!/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

   echo ERROR passing positions from LAMMPS to PLUMED

   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, $2-$6, $3-$7, $4-$8}'

fi

# CHECK PLUMED timestep

tstep=`grep timestep in.peptide-plumed | awk '{print $2}'`

tstep=`echo $tstep \* 0.001 \* 10 | bc -l`

nlines=`wc -l plmd_energy | awk '{print $1}'`

for ((i=3;i<$nlines;i++)); do

   told=`head -n $(($i-1)) plmd_energy | tail -n 1 | awk '{print $1}'`

   tnew=`head -n $i plmd_energy | tail -n 1 | awk '{print $1}'`

   tdiff=`echo \( $tnew - $told - $tstep \) \> 0 | bc -l`

   if [ $tdiff -gt 0 ] ; then

      echo ERROR passing timestep from LAMMPS to PLUMED

   fi

done

# Check PLUMED energy

tail -n +2 plmd_energy > plmd_energy2 

nlines=`paste lammps_energy plmd_energy2 | 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 potential energy from LAMMPS to PLUMED

   paste lammps_energy plmd_energy2 | tail -n +2 | awk '{if( $2<$4-0.0001 || $2>$4+0.0001 ) print $0}'

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

# Now run calculations to test virial

$LAMMPS < in.peptide-plumed-npt

$LAMMPS < in.peptide-plumed-npt2	

nlines=`paste plmd_volume_with_restraint plmd_volume_without_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 virial from PLUMED back to LAMMPS

fi

# Nothing from here works

# Now run calculations to check forces on energy

$LAMMPS < in.peptide-plumed-engforce-ref

$LAMMPS < in.peptide-plumed-eng-force-plumed

#!/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

#! FIELDS time dd

 0.000000 5.985554

 0.002000 6.002880

 0.004000 6.022015

 0.006000 6.029922

 0.008000 6.020103

 0.010000 5.996906

 0.012000 5.972734

 0.014000 5.960079

 0.016000 5.963714

 0.018000 5.978140

 0.020000 5.991813

 0.022000 5.995155

 0.024000 5.987021

 0.026000 5.975340

 0.028000 5.971456

 0.030000 5.981945

 0.032000 6.003550

 0.034000 6.024743

 0.036000 6.032990

 0.038000 6.022936

 0.040000 6.000131

 0.042000 5.977800

 0.044000 5.968692

 0.046000 5.977224

 0.048000 5.996934

 0.050000 6.014800

 0.052000 6.019586

 0.054000 6.008803

 0.056000 5.989809

 0.058000 5.974484

 0.060000 5.971140

 0.062000 5.979074

 0.064000 5.989379

 0.066000 5.991356

 0.068000 5.980176

 0.070000 5.960625

 0.072000 5.944401

 0.074000 5.942614

 0.076000 5.958402

 0.078000 5.984574

 0.080000 6.007964

 0.082000 6.017667

 0.084000 6.011795

 0.086000 5.998304

 0.088000 5.989405

 0.090000 5.993275

 0.092000 6.008545

 0.094000 6.025183

 0.096000 6.031186

 0.098000 6.020651

 0.100000 5.997952

 0.102000 5.975230

 0.104000 5.964757

 0.106000 5.971150

 0.108000 5.988568

 0.110000 6.004676

 0.112000 6.008731

 0.114000 5.998481

 0.116000 5.981406

 0.118000 5.969615

 0.120000 5.971827

 0.122000 5.987658

 0.124000 6.007888

 0.126000 6.020477

 0.128000 6.018377

 0.130000 6.004046

 0.132000 5.987682

 0.134000 5.980338

 0.136000 5.986534

 0.138000 6.001303

 0.140000 6.013589

 0.142000 6.013717

 0.144000 6.000028

 0.146000 5.980283

 0.148000 5.966836

 0.150000 5.968670

 0.152000 5.985459

 0.154000 6.007612

 0.156000 6.022374

 0.158000 6.022034

 0.160000 6.008851

 0.162000 5.993355

 0.164000 5.987212

 0.166000 5.995452

 0.168000 6.013111

 0.170000 6.028386

 0.172000 6.030387

 0.174000 6.016468

 0.176000 5.994191

 0.178000 5.976616

 0.180000 5.973983

 0.182000 5.987185

 0.184000 6.007275

 0.186000 6.021338

 0.188000 6.020837

 0.190000 6.006955

 0.192000 5.989433

 0.194000 5.979796

 0.196000 5.983601

 0.198000 5.996921

 0.200000 6.009310

 0.202000 6.011114

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

# 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