Update 4.FixTopoII/testing/test_serial/test6/run.lam.test.without,...

Update 4.FixTopoII/testing/test_serial/test6/run.lam.test.without, 4.FixTopoII/testing/test_serial/test6/run.lam.test.with, 4.FixTopoII/testing/test_serial/test6/README.md, 4.FixTopoII/testing/test_serial/test6/Eq.5_1.N1000.data, 4.FixTopoII/testing/test_serial/test6/bash.run.sh files

4.FixTopoII/testing/test_serial/test6/README.md

+# Test 6 in serial: The beads in the topoII section change type (from 1 to 2). But type 2 interacts in the same way as type 1.
+
+In this test we run two simulations:
+
+* 1.- Without introducing a topoII into the system (**run.lam.test.without**): In this case we use simply *fix nve + fix langevin* to run a Brownian Dynamics simulation. Only beads of type 1 are present into the system.
+* 2.- Introducing a topoII into the system (**run.lam.test.with**): In this case, besides the two previous fixes we use  *fix_topo2* with the following parameters:
+  * a) *n = 1*: the number of topoII regions.
+  * b) *l = 50 σ*: is the length of consecutive beads defining the topoII section. 
+  * c) *type = 2*: the beads in the section where the **topoII is jumping to** will be assigned this atom type.
+  * d) *seed:* to generate randomly the position of the topoII section
+  * e1) For the **static** keyword of this fix, in addition to (a)-(d) we need to set at which timestep the topoII will be introduced, this is *tin = 50*.
+  * e2) For the **randjump** and **jump2maxdens** keywords of this fix, in addition to (a)-(d) we need to set :
+    * The frequency at which the fix will be called (*freq=10*) 
+    * The probability that when the fix is called the topoII section changes position (*prob=1.0*). 
+    * The bead type that it will be temporally assigned to the section that the **topoII is jumping from** (*type=3*).
+
+Both simulations were run with the exact same seed for the Langevin fix to be able to compare the trajectory generated in the two cases. Also, both simulations start from the same initial configuration.
+
+**Note 1:** Since the *fix_topo2* changes the atom type of the particles in the section where it is jumping to (to type 2) and the section is jumping from (to type 3); then because type 2 and type 3 particles interact sterically in the same way as type 1 particles, we expect that neither velocities or forces are affecte by this change. If this is the case when running **bash.run.sh** at the end the programm will print success.
+
+**Note 2:** By default the lammps script with the topo2 (*run.lam.test.with*) sets the simulation with the keyword static. If you want to test a different keyword uncomment the lines related to fix_topo2 at the end of the file accordingly.

4.FixTopoII/testing/test_serial/test6/bash.run.sh

+#!/bin/bash
+
+#Directories
+cwd=$(pwd)
+
+#Executable
+LAMMPS=~/lmp_serial_29Oct2020_fixtopo2
+
+
+#Set variables
+initial=0
+last=10000
+dumpfreq=100
+check=0
+
+#Run system without fix topo2:
+  echo "Running system without TopoII"
+  cd ${cwd}
+  if [ -d "DUMPS_without_topo2" ]; then rm -Rf DUMPS_without_topo2; fi
+  
+  ${LAMMPS} -in run.lam.test.without
+  mv energy_total.txt DUMPS/
+  mv radius_of_gyration.txt DUMPS/
+  mv log.lammps DUMPS/log.lammps
+  mv DUMPS DUMPS_without_topo2
+
+
+#Run system with fix topo2:
+  echo "Running system with TopoII"
+  cd ${cwd} 
+  if [ -d "DUMPS_with_topo2" ]; then rm -Rf DUMPS_with_topo2; fi
+  
+  ${LAMMPS} -in run.lam.test.with
+  mv energy_total.txt DUMPS/
+  mv radius_of_gyration.txt DUMPS/
+  mv log.cite DUMPS/
+  mv log.lammps DUMPS/log.lammps1
+  mv DUMPS DUMPS_with_topo2
+
+
+#Compare results of the two simulations
+#create an empty file:
+if [ -f "results.txt" ] ; then
+    rm "results.txt"
+fi
+touch results.txt
+
+for t in $( seq 0 ${dumpfreq} ${last} )
+do
+  #cd ${cwd} 
+  
+  echo "timestep $t" >> results.txt
+
+
+  #Files containing: id type mol mass x y z ix iy iz vx vy vz fx fy fz
+  sim1="DUMPS_without_topo2/Run.R.${t}"
+  sim2="DUMPS_with_topo2/Run.R.${t}"
+  
+  #Remove the first 10 lines of each file
+   sed -i '1,9d' ${sim1}
+   sed -i '1,9d' ${sim2}
+
+  #Compare all the columns except ($2: type).
+  #FNR==NR is performed when reading the first file.
+  #{a[]; next}:       stores in a[] the columns of the first file and goes to the next line.
+  #($1,$3,...) in a:  is evaluated when looping through the second file. It checks if the current line is within the a[] array
+  #If all the columns (except column 2, related to the type) match between the two files, then print the first and second column of the second file and append it to results.txt
+  awk 'NR==FNR{a[$1,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12,$13,$14,$15,$16];next} ($1,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12,$13,$14,$15,$16) in a { print $1,$2;}' ${sim1} ${sim2} >> results.txt
+done
+
+#If all the columns match between the two files, then 11 lines (10 particles + 1 timestep) are printed per timestep into the file "results.txt".
+#number of particles
+np=1000
+#number of lines per timestep
+nl=$(( np + 1 ))
+#Total number of lines if all the columns matched for all timesteps: (11*101=1111)
+ntotal=$(( (np+1)*(((last-initial)/dumpfreq)+1) ))
+
+#Now count the number of lines in the file "results.txt"
+nlf=$(< "results.txt" wc -l)
+
+
+echo ${ntotal}
+echo ${nlf}
+
+if [ $ntotal -eq $nlf ]; then
+  echo "#####################"
+  echo "   TEST 6: SUCCESS   "
+  echo "#####################"
+fi
+
+if [ $ntotal -ne $nlf ]; then
+  echo "##################"
+  echo "   TEST 6: FAIL   "
+  echo "##################"
+fi

4.FixTopoII/testing/test_serial/test6/run.lam.test.with

+#################################
+####    PARAMETERS       ########
+#################################
+#temperature of the system
+variable Temp equal 1.0
+
+variable tdamp equal 1.0
+
+#Equilibration
+variable run1 equal 10000
+
+#Restart every these timesteps
+variable trestart equal 10000
+
+#Dump data every these timesteps
+variable dumpfreq equal 100
+
+#Dump trajectory every these timesteps
+variable dumptraj equal 100
+
+variable seedthermo equal 125485
+variable seedtopo   equal 79654
+
+################################
+####	DEFINTIONS 	    ########
+################################
+#Where the restart files of the run (non the equilibration) will be stored
+variable folder index DUMPS
+shell    mkdir ${folder}
+
+#The starting configuration is always the same. Also the name of the restart files
+variable rname   index Eq.5_1.N1000.data
+variable simname index Run.R
+
+units lj
+atom_style angle 
+boundary        p p p
+neighbor 1.2 bin
+neigh_modify every 1 delay 1 check yes 
+
+restart ${trestart} ${folder}/Restart.${simname}.
+read_data ${rname}
+
+
+####################################
+####	ANGLE	             #######
+####################################
+angle_style   cosine
+angle_coeff  1 5.0 #20 sigma for realistic DNA (1 bead = 2.5 nm)
+
+#######################################################
+####	PAIRS -- REPULSIVE + TOPO II 	#######
+#######################################################
+pair_style  hybrid lj/cut 1.12246 soft 1.12246
+pair_modify     shift yes
+pair_coeff      * * lj/cut 1.0 1.0 1.12246
+## NOW THERE IS TOPO2 -- TYPE 4 BEADS ARE SOFT AND CAN BE CROSSED
+pair_coeff      * 4 soft 2.0  1.12246
+pair_coeff      * 5 soft 20.0 1.12246
+
+###################################
+####  FENE	 ##################
+##################################
+bond_style   fene
+#It turns off the LJ part of the FENE-bond potential between 1-2 atoms. However, the pairwise LJ interaction is still on
+special_bonds fene#<=== I M P O R T A N T (new command).
+bond_coeff    1  30.0   1.6  1.0  1.0
+bond_coeff    2  30.0   1.6  1.0  1.0
+
+####################################
+####	THERMO AND INTEGRATION
+####################################
+thermo 1000
+thermo_style   custom   step  temp  epair vol cpu
+timestep 0.01
+
+####################################
+####	DUMPS	             #######
+####################################
+dump 2 all custom ${dumpfreq} ${folder}/${simname}.* id type mol mass x y z ix iy iz vx vy vz fx fy fz
+dump_modify 2 sort id
+
+##XYZ
+dump 4 all xyz ${dumptraj} ${folder}/traj.xyz
+dump_modify 4 element O N Ca P Au
+
+compute 	MSDs all msd
+compute		GYRs all gyration
+fix imprimir all print ${dumpfreq} "$(step)  $(c_GYRs) $(c_MSDs[4])" screen no append radius_of_gyration.txt
+
+#Compute total energies
+compute ekin all ke
+compute epot all pe
+variable ekin equal c_ekin
+variable epot equal c_epot
+variable etot equal c_ekin+c_epot
+fix energia all print ${dumpfreq} "$(step) ${ekin} ${epot} ${etot}" screen no append energy_total.txt
+
+####################################
+####    FIXES                #######
+####################################
+#fix 0 all topoII static        1 50 2 ${seedtopo} 100
+#fix 0 all topoII randjump      1 50 2 ${seedtopo} 100 1.0 3
+#fix 0 all topoII jump2maxdens  1 50 2 ${seedtopo} 100 1.0 3
+fix 0 all topoII jump2maxcurv  1 50 2 ${seedtopo} 100 1.0 3
+fix 1 all  nve
+fix 2 all  langevin   ${Temp} ${Temp} ${tdamp}  ${seedthermo}
+
+run ${run1}

4.FixTopoII/testing/test_serial/test6/run.lam.test.without

+#################################
+####    PARAMETERS       ########
+#################################
+#temperature of the system
+variable Temp equal 1.0
+
+variable tdamp equal 1.0
+
+#Equilibration
+variable run1 equal 10000
+
+#Restart every these timesteps
+variable trestart equal 10000
+
+#Dump data every these timesteps
+variable dumpfreq equal 100
+
+#Dump trajectory every these timesteps
+variable dumptraj equal 100
+
+variable seedthermo equal 125485
+
+################################
+####	DEFINTIONS 	    ########
+################################
+#Where the restart files of the run (non the equilibration) will be stored
+variable folder index DUMPS
+shell    mkdir ${folder}
+
+#The starting configuration is always the same. Also the name of the restart files
+variable rname   index Eq.5_1.N1000.data
+variable simname index Run.R
+
+units lj
+atom_style angle 
+boundary        p p p
+neighbor 1.2 bin
+neigh_modify every 1 delay 1 check yes 
+
+restart ${trestart} ${folder}/Restart.${simname}.
+read_data ${rname}
+
+
+####################################
+####	ANGLE	             #######
+####################################
+angle_style   cosine
+angle_coeff  1 5.0 #20 sigma for realistic DNA (1 bead = 2.5 nm)
+
+#######################################################
+####	PAIRS -- REPULSIVE + TOPO II 	#######
+#######################################################
+pair_style  hybrid lj/cut 1.12246 soft 1.12246
+pair_modify     shift yes
+pair_coeff      * * lj/cut 1.0 1.0 1.12246
+## NOW THERE IS TOPO2 -- TYPE 4 BEADS ARE SOFT AND CAN BE CROSSED
+pair_coeff      * 4 soft 2.0  1.12246
+pair_coeff      * 5 soft 20.0 1.12246
+
+###################################
+####  FENE	 ##################
+##################################
+bond_style   fene
+#It turns off the LJ part of the FENE-bond potential between 1-2 atoms. However, the pairwise LJ interaction is still on
+special_bonds fene#<=== I M P O R T A N T (new command).
+bond_coeff    1  30.0   1.6  1.0  1.0
+bond_coeff    2  30.0   1.6  1.0  1.0
+
+####################################
+####	THERMO AND INTEGRATION
+####################################
+thermo 1000
+thermo_style   custom   step  temp  epair vol cpu
+timestep 0.01
+
+####################################
+####	DUMPS	             #######
+####################################
+dump 2 all custom ${dumpfreq} ${folder}/${simname}.* id type mol mass x y z ix iy iz vx vy vz fx fy fz
+dump_modify 2 sort id
+
+##XYZ
+dump 4 all xyz ${dumptraj} ${folder}/traj.xyz
+dump_modify 4 element O N Ca P Au
+
+compute 	MSDs all msd
+compute		GYRs all gyration
+fix imprimir all print ${dumpfreq} "$(step)  $(c_GYRs) $(c_MSDs[4])" screen no append radius_of_gyration.txt
+
+#Compute total energies
+compute ekin all ke
+compute epot all pe
+variable ekin equal c_ekin
+variable epot equal c_epot
+variable etot equal c_ekin+c_epot
+fix energia all print ${dumpfreq} "$(step) ${ekin} ${epot} ${etot}" screen no append energy_total.txt
+
+####################################
+####    FIXES                #######
+####################################
+fix 1 all  nve
+fix 2 all  langevin   ${Temp} ${Temp} ${tdamp}  ${seedthermo}
+
+run ${run1}