[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[ESPResSo] Re: non-bonded interactions not showing up in the espresso po
|
From: |
Christopher Jesudason |
|
Subject: |
[ESPResSo] Re: non-bonded interactions not showing up in the espresso polymer script for analyze energy command |
|
Date: |
Fri, 24 Aug 2007 06:33:42 -0700 (PDT) |
Dear friends,
I have a very simple single polymer electrolyte with
species type 0 (150 of them ) all charged -1
interacting with counterions type 1 all charged +1
[with harmonic potential type o and bending potential
type 1 for this polymer chain]. When I used the
analyze command to get the nonbonded interaction
energies, where a repulsive LJ potential is swithced
on for the 0 0 , 1 0 and 1 1 ions only the (0 0)
energies show up , but the energy for the (1 0) and (0
0) non-bonded interactions of LJ type both show zero
in both cases, which is impossible for the case of the
(0 1 ) interactions because they have opposite
charges.
I looked into the manual and there is something in
yellow on page 59 (top lhs) with comments that cannot
be read because it was placed outside the text
margin, at least according to the viewer used here.
the snippet of the output is :
all energy { energy -155.547420243 } { kinetic
471.398271936 } { 0 HARMONIC 113.941433298 } { 1 angle
0.0 } { 0 0 nonbonded 54.5179976882 } { 0 1
nonbonded 0.0 } { 1 1 nonbonded 0.0 } { coulomb
-795.405123166 0.0 -795.405123166 }
all total -155.547420243
linindex01 -155.547420243
linindex11 471.398271936
linindex22 113.941433298
kinetic 471.398271936
nonbonded 00 54.5179976882
nonbonded 10 0.0
nonbonded 11 0.0
113.941433298
bonded 113.941433298
coulomb -795.405123166
Enclosed the code . I would be very pleased if you
could give me an indication of why this is the case.
Has there been an ommision of the energies or is the
total LJ interaction energy put into the 0 0
interaction index in analyse energy command?
Please clarify. And thank you for your time. I hope I
will be able to complete my research project with
Espressso code. Could it be that the myconfig file
needs to be activated?
Chis
Chris
--- Olaf Lenz <address@hidden> wrote:
>
>
____________________________________________________________________________________
Boardwalk for $500? In 2007? Ha! Play Monopoly Here and Now (it's updated for
today's economy) at Yahoo! Games.
http://get.games.yahoo.com/proddesc?gamekey=monopolyherenow #!/bin/sh
# tricking... the line after a these comments are interpreted as standard shell
script \
exec $ESPRESSO_SOURCE/Espresso $0 $*
puts "[code_info]"
set syst1 [clock seconds]
#exit
# random numbers intialization
proc RandomInit { seed } {
global randomSeed
set randomSeed $seed
}
proc Random {} {
global randomSeed
set randomSeed [expr ($randomSeed*9301 + 49297) % 233280]
return [expr $randomSeed/double(233280)]
}
proc RandomRange { range } {
expr int([Random]*$range)
}
# Add configuration in XMOL format to the trajectory file "f_xmol"
proc confxmol { f_xmol } {
global box_l n_part chem_typ
puts $f_xmol [expr int($n_part)]
puts $f_xmol " after [setmd time] i.u., BOX: $box_l $box_l $box_l"
for {set j 0} {$j < $n_part} {incr j} {
puts $f_xmol " $chem_typ($j) [part $j print pos] "
}
}
;# left out the gromacs file output
# Put molecules into the central box if they go out
# Array moladr($n_mol) must be defined (= first atom of each molecule)
#
proc ToPBC {} {
global box_l n_mol moladr
# loop over molecules
for {set i 0} {$i<$n_mol} {incr i} {
# compute center-of-mass for each molecule
set xcom 0
set ycom 0
set zcom 0
set ii [expr $i + 1]
## puts "$i $ii $moladr($i) $moladr($ii)"
for {set j $moladr($i)} {$j < $moladr($ii)} { incr j } {
set xcom [expr $xcom + [lindex [part $j print pos] 0 ] ]
set ycom [expr $ycom + [lindex [part $j print pos] 1 ] ]
set zcom [expr $zcom + [lindex [part $j print pos] 2 ] ]
}
set xcom [ expr $xcom / ($moladr($ii)-$moladr($i))]
set ycom [ expr $ycom / ($moladr($ii)-$moladr($i))]
set zcom [ expr $zcom / ($moladr($ii)-$moladr($i))]
# find eventual shift
if { $xcom < 0. } then {
set xshift $box_l
} elseif { $xcom > $box_l} {
set xshift [expr -1*$box_l]
} else {
set xshift 0.
}
if { $ycom < 0. } then {
set yshift $box_l
} elseif { $ycom > $box_l } {
set yshift [expr -1*$box_l]
} else {
set yshift 0.
}
if { $zcom < 0. } then {
set zshift $box_l
} elseif { $zcom > $box_l } {
set zshift [expr -1*$box_l]
} else {
set zshift 0
}
# move the whole molecule so that its COM is inside the main cell
for {set j $moladr($i)} {$j < $moladr($ii)} { incr j } {
set xnew [expr [lindex [part $j print pos] 0] + $xshift]
set ynew [expr [lindex [part $j print pos] 1] + $yshift]
set znew [expr [lindex [part $j print pos] 2] + $zshift]
part $j pos $xnew $ynew $znew
}
}
}
#puts "goodbye"
#exit
#############################################################
#
proc ToCENTER { } {
global box_length n_part ifm ilm
# assumes the same mass for each gorup of particles in a molecule etc..
set xcom 0
set ycom 0
set zcom 0
# set ii [expr $i + 1]
# loop over molecules
for {set j $ifm } {$j <= $ilm} {incr j} {
# compute center-of-mass for each molecule
## puts "$i $ii $moladr($i) $moladr($ii)"
set xcom [expr $xcom + [lindex [part $j print pos] 0 ] ]
set ycom [expr $ycom + [lindex [part $j print pos] 1 ] ]
set zcom [expr $zcom + [lindex [part $j print pos] 2 ] ]
}
set xcom [ expr $xcom / ( $ilm-$ifm +1)]
set ycom [ expr $ycom /( $ilm-$ifm +1) ]
set zcom [ expr $zcom /( $ilm-$ifm +1) ]
###puts "shifting DNA COM to the center"
for {set j 0} {$j < $n_part} { incr j } {
set xnew [expr [lindex [part $j print pos] 0] - $xcom +
0.5*$box_length]
set ynew [expr [lindex [part $j print pos] 1] - $ycom +
0.5*$box_length]
set znew [expr [lindex [part $j print pos] 2] - $zcom +
0.5*$box_length]
part $j pos $xnew $ynew $znew
}
;##puts "Tcenter $j "
}
#puts "goodbye"
#exit
;# we place the snake integration parameters here
# 1.6 Integration parameters
set time_step 0.01
set skin 0.4
set gamma 0.5
set temp 1.00
# here placed the thermostat command for NVT type
####### integrate set nvt
thermostat langevin $temp $gamma
setmd time_step $time_step
setmd skin $skin
# integration type is also mentioned
## warm-up integration
# 1.6.1 warmup integration (with capped LJ potential) until particles are at
least $min_dist apart
# warm steps in actual run is 100 and warm loops at 200
set warm_steps 10
set warm_loops 200
set warm_cap1 20
set warm_incr 2
# 1.6.2 integration (with full LJ potential) for $int_time
# actual cgj set int_steps 100 for full runs use these parameters
#actual cgj set int_times 10000 for full runs use these parameters below
only for tests.
# below tests only
set int_steps 100
set int_times 40000
### dmp values
#no of times for printout (200 is the norm)
set dmp 100
puts " number of dumps is : $dmp "
##require to set dmp interval
set dmpintv [expr $int_times/$dmp]
puts "dump interval is $dmpintv"
set dmpfac [expr 1.00/double($dmpintv)]
;# set up the box
# 2. Setup System
# 2.1 setting the box and all interactiion parameters:
puts " "
puts "==================================================="
puts "= p150.tcl polymer ="
puts "==================================================="
puts " "
# 1.2 System parameters
# Units: length = A
# k = 1
# m = 1 i.u.
# charge = -e
# T = 1
# E in kT units, that is E=1 is 2.4 kJ/M at 300k
# 1.3 Particles
# Num of NCPs counter ion particles
set n_C 150
# Charge and radius of Core
set q_C 1.
set rO 0.0
# Num of DNA CG units
set n_D 150
# Charge and radius of CG DNA blobs (blob = 6 bp)
set Q_D -1.
set R_D 0.0
# bond length
set R_link 1.00
# calculated from neutral + salt concentration?
#set nNa 2240
#set qNa 1
#set nCl 100
#set qCl -1
# set masses (masC- NCP, masD - DNA blobs)
set masC 1.
set masD 1.
set mass 1.
#set masNa 1.
#set masCl 1.
# Box size
set box_length 300.0
set box_l 300.
set shield 1.
# Interaction parameters here:
set n_part_types 2
# repulsive Lennard Jones : all cut-off at power(2,1/6)*sigma : shift
1/4*epsilon : sigma radius
# for central units
set ljC_eps 1.
set ljC_sig 1.
set ljC_cut [expr 1.12246*$ljC_sig]
set ljC_shift [expr 0.25*$ljC_eps]
set ljC_off 0.
### $rO before this was ro
# for DNA
set ljD_eps 1.
set ljD_sig 1.
set ljD_cut [expr 1.12246*$ljD_sig]
set ljD_shift [expr 0.25*$ljD_eps]
set ljD_off 0.
### bwfor this was rd $R_D
# for ions
# some parameters
#set ljI_eps 1.
#set ljI_sig 2.5
#set ljI_cut [expr 1.12246*$ljI_sig]
#set ljI_shift [expr 0.25*$ljI_eps]
#set ljI_off 0.0
# harmonic stiffer bonds
set harm_k 100.00
set harm_r $R_link
set harm_k [expr 100.00*$temp]
# angle bending (DNA persistence length)
set bend_k 0.
# leafyoung
# 1.5 Electrostatics
# Bjerrum length corresponding to water dielectric
# Tuning
set p3m_bjerrum 1.00
set p3m_accuracy 9.9722e-07
# After tuning
set p3m_r_cut 3.18484e-01
set p3m_mesh 16
set p3m_cao 5
set p3m_alpha 8.23081e+00
# file parameters and labels
# 1. Parameters
# 1.1 System identification:
set name "p150"
set ident "v1"
# read restart?
set read_conf no
# 1.2 System parameters
# Files for observables:
set f_xmol [open polymer.xmol "w"]
set f_gro [open spid.gro "w"]
set f_conn [open conn.tcl "w"]
set rdf_obs [open rdf.dat "w"]
## SCRATCH FILE
#####set rdf_scr [open rdf.dat "w"]
set f_scr [open "$name$ident.dat" "w"]
;# 1.7 restart point
set restart poly.dmp
##set rrestart poly.chk
# observables a file for observables
set f_obs [open "$name$ident.obs" "w"]
set fwarm_obs [open "w$name$ident.obs" "w"]
;# now puts the number of units you desire
set l_poly $n_D
set n_ci $n_C
set n_part [expr $l_poly+$n_ci]
### setmd n_part $n_part
puts "Simulate PE Solution N=$l_poly at density due to chain and counterions"
puts "Simulation box: $box_length"
;# 1.7 restart point
####set restart "poly.dmp"
#######set rrestart "poly.chk"
;#puts "haha[ expr $box_length]"
;#exit
setmd box_l $box_length $box_length $box_length
#puts " [setmd box_l]"
#exit
#############################################################
#
# We have to introduce the electrostatic interactions #
# later for technical reasons. #
#############################################################
# "connection file" is used by gOpenMol to draw a nicer picture
# opening file for gOpenMol we start the f_conn file at this point
puts $f_conn "atom cpk * * *"
# i - particle counter
set i 0
# im - molecular counter
set im 0
set im1 0
# molecular addresses
set moladr(0) 0
# require to set molecular configuration and addresses
puts $f_conn "atom scale cpk 8. * * $i"
# 2.2.1 Define Bond type 0 (between Tails), 1 (between O-Tail)
inter 0 harmonic $harm_k $harm_r
##exit
inter 1 angle $bend_k
;# cgj orig inter 1 angle $bend_k
;############## inter 1 BOND_ANGLE_HARMONIC $bend_k
# 2.3 DNA
puts "Generating DNA of length $n_D charge $Q_D"
set nd1 [expr $n_D + 1]
polymer 1 $l_poly $R_link start $i charge $Q_D distance 1 types 0 0
# loop over DNA units
for {set ii 0} {$ii < $n_D} {incr ii} {
part $i type 0 mass $masD q $Q_D
set chem_typ($i) "C"
set i1 [expr $i + 1]
if {$i1 == $n_D } break
puts $f_conn "atom scale cpk 8. * * $i1"
set ic [expr $i - 1]
if { $ii == 0} then {
puts "start DNA from random position"
} else {
part $ic bond 0 $i
puts $f_conn "edit bond create * * $ic * * $i"
if { [expr $ii + 1] < $n_D} {
part $ic bond 1 $i $i1
}
}
# for loop close after incri
incr i
}
incr im
set moladr($im) [expr $moladr($im1)+$n_D]
## moladr(1) is the first particle of the counter ions Na
# f_conn defined here
#######polymer 1 $l_poly 1.0 start 0 charge 1.0 types 0 0 FENE 0
set j 0
counterions $n_ci start $l_poly charge $q_C type 1
for {set j $l_poly} {$j <[expr $n_part ]} {incr j} {
# might want to add mass here
part $j type 1 q $q_C mass $masC
puts $f_conn "atom scale cpk 8. * * $j"
set chem_typ($j) "O"
# must also increment moladr
set moladr($im) $j
incr im
}
puts " im $im "
# nb im always points one ahead
puts [part 299]
####puts [part 300]
puts [part 99]
#exit
#puts "[inter coulomb 1.0 p3m tune accuracy 0.001 mesh 8]"
#Defining electrostatics just before warm-up
##puts "first [inter coulomb 1.00 p3m tunev2 accuracy 1.e-6 mesh 32]"
# for this 150 system
#resulting parameters:
#16 5 4.18654e-01 6.80699e+00 9.99703e-07 9
#####################
###inter coulomb 1.00 p3m tune accuracy 1.e-5 mesh 32
puts "tuning here "
#puts "[inter coulomb 1.0 p3m tune accuracy 1.e-6 mesh 32]"
#####puts "second [inter coulomb 1.00 p3m tunev2 accuracy 1.e-6 r_cut 0 mesh
0 cao 0]"
inter coulomb $p3m_bjerrum p3m $p3m_r_cut $p3m_mesh $p3m_cao $p3m_alpha
$p3m_accuracy
# just before the box setup
# 2.2 setting up interactions
# ------------------------
;# inter 1 angle $bend_k
#puts " here"
inter 0 0 lennard-jones $ljD_eps $ljD_sig $ljD_cut $ljD_shift [expr 2*$ljD_off]
##inter 2 2 lennard-jones $ljC_eps $ljC_sig $ljC_cut $ljC_shift [expr
2*$ljC_off]
inter 0 1 lennard-jones $ljD_eps $ljD_sig $ljD_cut $ljD_shift [expr $ljC_off +
$ljD_off ]
inter 1 1 lennard-jones $ljC_eps $ljC_sig $ljC_cut $ljC_shift [expr
2*$ljC_off]
####inter 1 1 lennard-jones $ljD_eps $ljD_sig $ljD_cut $ljD_shift [expr
2*$ljD_off]
##exit
;## here used the older warm up criterion
# ifm is the first subscript for molecule
# ilm is the second subscript for molecule
# Set LJ cup
set ifm 0
set ilm [ expr $l_poly -1]
set n_mol [expr $im -1]
# due to extra factor please modify this accordingly for other molecular types
#puts " nmol number of molecular species $n_mol"
##puts " npartalready defined $n_part"
# due to extra factor please modify this accordingly for other molecular types
ToCENTER
ToPBC
puts "here in center"
###exit
### need to do a ToPBC here
# Dump file with start-up in Espresso format:
###### polyBlockWrite "$name$ident.0000" all
confxmol $f_xmol
close $f_conn
# read restart coordinates
if { $read_conf == yes } then {
puts "trying to read restart file $restart"
#####checkpoint_read $rrestart
set restart [open $restart "r"]
blockfile $restart read auto
blockfile $restart read particles
close $restart
puts "Old configuration restarted successfully!"
set warm_steps 0
set warm_loops 0
;######################################set f_trxmol [open snak.xmol "w"]
} else {
setmd time 0
}
puts " Ready to start"
puts "traj file $f_xmol"
set cap $warm_cap1
inter ljforcecap $cap
# Warmup Integration Loop
set i 0
while { $i < $warm_loops } {
integrate $warm_steps
# Warmup criterion
set act_min_dist [analyze mindist]
set act_min_dist_c [analyze mindist 1 1]
set act_energy [analyze energy]
set bond_ener [expr [analyze energy bonded 0]+[analyze energy bonded 1]]
set el_energy [analyze energy coulomb]
puts "t=[setmd time] LJcap=$cap rm= $act_min_dist_c \
E= [lindex $act_energy 0 1] $bond_ener $el_energy\r"
flush stdout
# write observables
puts $fwarm_obs "{ time [setmd time] } $act_energy"
# Increase LJ cap
set cap [expr $cap+$warm_incr]
inter ljforcecap $cap
incr i
# put configuration in XMOL trajectory; need to check on this format first.
## ToPBC
;# ,,, this needs to be modified please check the configuration
#
##confxmol $f_xmol
## ;;need to be changed
#xtc_write
}
ToPBC
confxmol $f_xmol
integrate $int_steps
puts "Warmup finished."
puts "all energy [analyze energy]"
puts "all total [analyze energy total]"
set act_energy [analyze energy]
puts "linindex01 [lindex $act_energy 0 1] "
puts "linindex11 [lindex $act_energy 1 1] "
puts "linindex22 [lindex $act_energy 2 2] "
puts " kinetic [analyze energy kinetic ]"
puts " nonbonded 00 [analyze energy nonbonded 0 0]"
puts " nonbonded 10 [analyze energy nonbonded 1 0]"
puts " nonbonded 11 [analyze energy nonbonded 1 1]"
puts "[expr [analyze energy bonded 0]+[analyze energy bonded 1]]"
puts " bonded [analyze energy bonded 0]"
puts "coulomb [analyze energy coulomb]"
exit
;# turn off the ljforcecap, which is done by setting the
# force cap value to zero:
inter ljforcecap 0
puts "setmd particles [setmd n_part]"
# Begin production run here ;;;
## We also put some scratch files here to see how things can happen
## Initialize rdf files
set rdfbin 100
set cnt 0
for {set i 0} {$i < $rdfbin } {incr i} { lappend avg_rdf12 0}
for {set i 0} {$i < $rdfbin } {incr i} { lappend avg_rdf22 0}
for {set i 0} {$i < $rdfbin } {incr i} { lappend avg_rdf11 0}
set rdf11 [analyze rdf 1 1 0.0 [expr $box_l/2] $rdfbin]
set rdf12 [analyze rdf 1 2 0.0 [expr $box_l/2] $rdfbin]
set rdf22 [analyze rdf 2 2 0.0 [expr $box_l/2] $rdfbin]
set rlist11 ""
set rdflist11 ""
set rlist22 ""
set rdflist22 ""
set rlist12 ""
set rdflist12 ""
## setting the rlist that need not be repeated
foreach value [lindex $rdf11 1] {
lappend rlist11 [lindex $value 0]
}
;#puts "avg_rdf=$avg_rdf"
foreach value [lindex $rdf12 1] {
lappend rlist12 [lindex $value 0]
}
;#puts "avg_rdf=$avg_rdf"
foreach value [lindex $rdf22 1] {
lappend rlist22 [lindex $value 0]
}
;#puts "avg_rdf=$avg_rdf"
puts "\nStart integration: run $int_times times $int_steps steps"
set i 0
set cnt 0
set dmpval 0
### setting initial values here
set s_etot 0.
set s_act_energy01 0.
set s_bond_ener 0.
set s_el_energy 0.
set s_tempr 0.
set s_rg 0.
set s_Pres 0.
set s_re0 0.
set s_re2 0.
set s_tempr2 0.
set s_kinetic 0.
set s_potl 0.
######################################################
while { $i < $int_times } {
integrate $int_steps
# this not calculated now
# set act_min_dist [analyze mindist]
# set act_min_dist_c [analyze mindist 1 1]
# radius of gyration here
analyze set chains 0 1 $l_poly
### this portion deals with the radial distribution function
###########################################################################################
set rdf11 [analyze rdf 1 1 0.0 [expr $box_l/4.] $rdfbin]
set rdf12 [analyze rdf 1 2 0.0 [expr $box_l/4.] $rdfbin]
set rdf22 [analyze rdf 2 2 0.0 [expr $box_l/4.] $rdfbin]
## must zero each file list
##set rlist11 ""
set rdflist11 ""
###### set rlist22 ""
set rdflist22 ""
####### set rlist12 ""
set rdflist12 ""
#######################################
foreach value [lindex $rdf11 1] {
lappend rdflist11 [lindex $value 1] }
set avg_rdf11 [vecadd $avg_rdf11 $rdflist11]
;#puts "current avg_rdf AND CNT= $cnt avg rdf = $avg_rdf "
;# oo interactions
foreach value [lindex $rdf12 1] {
lappend rdflist12 [lindex $value 1] }
set avg_rdf12 [vecadd $avg_rdf12 $rdflist12]
;# 11 interactions
foreach value [lindex $rdf22 1] {
lappend rdflist22 [lindex $value 1] }
set avg_rdf22 [vecadd $avg_rdf22 $rdflist22]
;#end of particle reading block puts "again avg_rdf=$avg_rdf"
;# puts "cnt= $cnt +1 to counter number of times sampling taken "
set cnt [expr $cnt + 1]
#################################################################################################
set rg [lindex [analyze rg] 0]
set re [analyze re]
# please check this carefully we need to put this also in a file for
observation
## puts $obs "[setmd time] $rg"
## further tests because of ambiguity of description of end to end distance
### set xdiff [expr [lindex [part 0 print pos] 0 ]-[lindex [part 299 print
pos] 0 ] ]
### set ydiff [expr [lindex [part 0 print pos] 1 ]-[lindex [part 299 print
pos] 1 ] ]
###set zdiff [expr [lindex [part 0 print pos] 2 ]-[lindex [part 299 print
pos] 2 ] ]
####set e2 [expr pow($xdiff,2)+pow($ydiff,2)+ pow($zdiff,2) ]
###### set eabs [expr sqrt($e2)]
set Pres [analyze pressure total]
set act_energy [analyze energy]
set etot [analyze energy total]
set bond_ener [expr [analyze energy bonded 0]+[analyze energy bonded 1]]
set el_energy [analyze energy coulomb]
set kin_energy [analyze energy kinetic]
set tempr [expr [lindex $act_energy 1 1]/(1.5*$n_part)]
set tempr2 [expr $kin_energy/(1.5*[setmd n_part])]
puts "t=[setmd time] E= [lindex $act_energy 0 1] Eb= $bond_ener Eel=
$el_energy T= $tempr Rg=$rg"
set s_etot [expr $s_etot + $etot ]
set s_act_energy01 [expr $s_act_energy01 + [lindex $act_energy 0 1] ]
set s_bond_ener [expr $s_bond_ener + $bond_ener ]
set s_el_energy [expr $s_el_energy + $el_energy ]
set s_tempr [expr $s_tempr + $tempr ]
set s_rg [expr $s_rg + $rg ]
set s_Pres [expr $s_Pres + $Pres ]
set s_re0 [expr $s_re0 +[lindex $re 0] ]
set s_re2 [expr $s_re2 + [lindex $re 2] ]
set s_tempr2 [expr $s_tempr2 + $tempr2 ]
;###puts "setmd particles [setmd n_part]"
# write observables
#### puts $f_obs "{ time [setmd time] } $act_energy $rg "
# puts $f_obs " [setmd time] [lindex $act_energy 0 1] $bond_ener
$el_energy $tempr $rg"
## puts $f_scr " [setmd time] $etot $Pres [lindex $re 0] [lindex $re 2]
$tempr2"
####### puts $f_scr " hand computed $e2 , $eabs "
if {[expr ($i+1)%$dmpintv]==0 } {
set dmpval [expr $dmpval + 1]
puts $f_obs " $dmpval [expr $dmpfac*$s_act_energy01] [expr
$dmpfac*$s_bond_ener] [expr $dmpfac*$s_el_energy]\
[expr $dmpfac*$s_tempr] [expr $dmpfac*$s_rg]"
puts $f_scr " $dmpval [expr $dmpfac*$s_etot] [expr $dmpfac*$s_Pres] [expr
$dmpfac*$s_re0]\
[expr $dmpfac*$s_re2] [expr $dmpfac*$s_tempr2]"
set s_etot 0.
set s_act_energy01 0.
set s_bond_ener 0.
set s_el_energy 0.
set s_tempr 0.
set s_rg 0.
set s_Pres 0.
set s_re0 0.
set s_re2 0.
set s_tempr2 0.
}
incr i
# put configuration in XMOL trajectory
ToPBC
confxmol $f_xmol
# at this place after incr do we do the dmp values
}
integrate $int_steps
#######invalidate_system
set restart [open $restart "w"]
blockfile $restart write variable time
blockfile $restart write particles
close $restart
puts "restart point written"
###### checkpoint_set $restart
###### invalidate_system
####blockfile $restart write
###polyBlockWriteAll "$restart" all
# I set the checkpoint outside the loops since it is anticipated that at most a
continuity is required
### final averaging for the rdf outside the loop
set avg_rdf11 [vecscale [expr 1.0/$cnt] $avg_rdf11]
set avg_rdf12 [vecscale [expr 1.0/$cnt] $avg_rdf12]
set avg_rdf22 [vecscale [expr 1.0/$cnt] $avg_rdf22]
puts $rdf_obs "\# r11 r12 r22 rdf11 rdf12 rdf22 "
foreach r11 $rlist11 r12 $rlist12 r22 $rlist22 rdf1 $avg_rdf11 rdf12
$avg_rdf12 \
rdf22 $avg_rdf22 { puts $rdf_obs "$r11 $r12 $r22 $rdf1 $rdf12 $rdf22" }
#temperory to test the look of the full rdf which was never properly described
; then need to
### average and printout the three rdf's before running them
### integrate $int_steps
####set rdft [analyze rdf 1 2 0.0 [expr $box_l/2] $rdfbin]
######puts $rdf_obs "$rdft "
#xtc_close
close $f_obs
close $fwarm_obs
close $f_xmol
close $f_gro
close $rdf_obs
close $f_scr
########close $restart
set syst2 [clock seconds]
set dif [expr $syst2- $syst1]
puts "The elapsed time is $dif seconds "
;##########close $f_conn
exit
#############################################################