Commit b1d87791 authored by Jony Castagna's avatar Jony Castagna

Merge branch 'Thermodynamic_Force_Calc' into 'master'

Thermodynamic force calc

See merge request e-cam/E-CAM-Library!74
parents 2766bed5 3f4869fd
#!/bin/bash
source /<our path here>/gromacs/bin/GMXRC
j=$1
l=$2
k=$3
m=$4
n=$5
rmin= "Start of interpolation zone"
rmax= "End of interpolation zone"
rbox= "total length of the box"
rref= "Reference coordinate used in the AdResS input"
lc= "how fine the grid for the density should be"
prefac= "either select the prefactor from the formula or define a value here"
rstep=$(echo "$rbox/($lc-1)" | bc -l)
r_cut=$rmin
rm_c=$rmax
echo $r_cut, $rm_c
prep_sys=$j
case $prep_sys in
[0]*)
if [ ! -f F_th_step0 ]; then
seq -f%g 0.000 $rstep $rbox | sed 's/,/./' | awk '{print $1, "0.0", "0.0"}'> F_th_step0
fi
cp F_th_step0 tabletf_WCG.xvg
gmx mdrun -v
;;
[1]*)
if [ -f tabletf_CAT.xvg ]; then
awk 'BEGIN{OFS="\t"}(NR>9){print $1, $2, $3}' tabletf_CAT.xvg > F_th_step0
fi
if [ -f tabletf_ANN.xvg ]; then
awk 'BEGIN{OFS="\t"}(NR>9){print $1, $2, $3}' tabletf_ANN.xvg > F_th_step0
fi
;;
[2]*)
if [ -f tabletf_WCG.xvg ]; then
awk 'BEGIN{OFS="\t"}(NR>9){print $1, $2, $3}' tabletf_WCG.xvg > F_th_step0
fi
;;
esac
#rn=$(tail -n 1 F_th_step0 | awk '{printf("%f\n",$1)}')
rn1=$(echo $rn+$rstep | bc -l)
echo $rn $rn1
dens_sys=$n
dens_chk=$n
case $dens_sys in
[0]*)
csg_density --axis r --rmax 9.6 --ref [5.5,5.5,5.5] --trj ../7000_mmin/traj_comp.xtc --top ../7000_mmin/topol.tpr --out test.dens.comp
csg_density --axis r --rmax 9.6 --ref [5.5,5.5,5.5] --cg "MIL.cg.xml;cl.cg.xml" --trj traj_comp.xtc --top topol.tpr --molname Cl --out Cl.dens.out
;;
[1]*)
source ~/votca_5.1/bin/VOTCARC.bash
csg_density --axis r --rmax 10. --ref [4.89439,4.89439,4.89439] --trj traj_comp.xtc --top topol.tpr --out SOL.dens.dat
sed 's/nan/0.0/' SOL.dens.dat > calc_dens
seq -f%g 0 $rstep $rbox | sed 's/,/./' | awk 'BEGIN{OFS="\t"}{print $1, "575.0"}' > ref_dens
;;
[2]*)
echo "2" > t_c; gmx density -d X -sl $lc -f traj_comp.xtc -o SOL.dens.xvg < t_c
awk 'BEGIN{OFS="\t"}(NR>24){print $1, $2}' SOL.dens.xvg > $k
seq -f%g 0 $rstep $rbox | sed 's/,/./' | awk 'BEGIN{OFS="\t"}{print $1, "1000.0"}' > $l
#folding and symmetrizing density for smoothing:
# give asbolute values back: z=${res/#-/}; echo $z
awk '{d=$1-'$rref'; print ((d>0)?d:-d), $2}' $l > ref_dens_t
awk '{d=$1-'$rref'; print ((d>0)?d:-d), $2}' $k > calc_dens_t
rm x1 x2 xaa xab ref_dens calc_dens
#lc1=$(wc -l ref_dens_t | awk '{print $1/2}')
#echo $lc1
split -l $(( ($lc-1)/2 )) ref_dens_t
awk '{printf "%10f %10f\n", $1, $2}' xaa | sort -g > x1
awk '{printf "%10f %10f\n", $1, $2}' xab | sort -g > x2
paste x1 x2 | awk '{print $1,($2+$4)/2.0}' > ref_dens
#awk '{print $1,$2}' x2 > ref_dens
split -l $(( ($lc-1)/2 )) calc_dens_t
awk '{printf "%10f %10f\n", $1, $2}' xaa | sort -g > x1
awk '{printf "%10f %10f\n", $1, $2}' xab | sort -g > x2
paste x1 x2 | awk '{print $1,($2+$4)/2.0}' > calc_dens
;;
esac
seq -f%g 0.0 $rstep $rbox | sed 's/,/./' > xnew_ntrpl
paste calc_dens ref_dens | awk '{if($1<('$r_cut')) print $1,($2-$4)}' > dens_mix_at
paste calc_dens ref_dens | awk '{if(($1>=('$r_cut')) && ($1<('$rm_c'))) print $1,$2,$4}' > dens_mix_hy
paste calc_dens ref_dens | awk '{if($1>=('$rm_c')) print $1,($2-$4)}' > dens_mix_cg
./smooth_dens.sh $prefac $rstep 400.
cp dens_smooth d_s
cp pot_smooth p_s
echo "#","manually gen. Thermodyn. Force approx." > $m.xvg
echo "#","Parameter:">> $m.xvg
echo "#","start hy-region:", $rmin >> $m.xvg
echo "#","end of hy-region:",$rmax >> $m.xvg
echo "#","start tf: from xsplit" >> $m.xvg
echo "#","start hy-region:", $rmin >> $m.xvg
echo "#","end of hy-region:",$rmax >> $m.xvg
echo "#" >> $m.xvg
echo "#" >> $m.xvg
fmax=$(head -n 1 d_s | awk '{print $2}')
fmin=$(tail -n 1 d_s | awk '{print $2}')
pmax=$(head -n 1 p_s | awk '{print $2}')
pmin=$(tail -n 1 p_s | awk '{print $2}')
r_hy_1=$(head -n 1 d_s | awk '{print $1-'$rstep'}')
r_hy_2=$(tail -n 1 d_s | awk '{print $1}')
rbox1=$(echo $rbox-$rstep | bc -l)
seq -f%g 0.000 $rstep $r_hy_1 | sed 's/,/./' | awk '{print $1, 0.0}' > d0
awk '{print $1, $2}' d_s > d1
seq -f%g $r_hy_2 $rstep $rbox1 | sed 's/,/./' | awk '{print $1, 0.0}' > d2
seq -f%g 0.000 $rstep $r_hy_1 | sed 's/,/./' | awk '{print $1, '$pmax'}' > p0
awk '{print $1, $2}' p_s > p1
seq -f%g $r_hy_2 $rstep $rbox1 | sed 's/,/./' | awk '{print $1, '$pmin'}' > p2
cat d0 d1 d2 > d_m
cat p0 p1 p2 > p_m
paste p_m d_m | awk '{print $1,$2,$4}' > dens_mix
paste F_th_step0 dens_mix | awk 'BEGIN{OFS="\t"}{printf("%e %e %e\n",$1,($2-$5),($3-$6)) }' >> $m.xvg
cp $m.xvg tabletf_WCG_$m.xvg
cp $m.xvg tabletf_WCG.xvg
#clean-up:
gmx mdrun -v
dens_chk=$n
case $dens_chk in
[0]*)
csg_density --axis r --rmax 9.6 --ref [5.5,5.5,5.5] --trj ../7000_mmin/traj_comp.xtc --top ../7000_mmin/topol.tpr --out test.dens.comp
;;
[1]*)
source ~/votca_5.1/bin/VOTCARC.bash
csg_density --axis r --rmax 10. --ref [4.89439,4.89439,4.89439] --trj traj_comp.xtc --top topol.tpr --out SOL.dens.xvg
sed 's/nan/0.0/' SOL.dens.xvg > SOL.dens_$m.xvg
;;
[2]*)
echo "2" > t_c; gmx density -d X -f traj_comp.xtc -o SOL.dens_$m.xvg < t_c
;;
esac
# CLEAN UP
rm p? d? x? x?? \#*
rm SOL.dpot.*
rm p_m d_m
rm d_s d_s_t
rm t_c
rm s?.xvg s??.xvg
rm ref_dens
rm ref_dens_t calc_dens*
rm s.d.o
rm p_s
rm ref_t
rm calc_t
rm dens.ref
:orphan:
.. In ReStructured Text (ReST) indentation and spacing are very important (it is how ReST knows what to do with your
document). For ReST to understand what you intend and to render it correctly please to keep the structure of this
template. Make sure that any time you use ReST syntax (such as for ".. sidebar::" below), it needs to be preceded
and followed by white space (if you see warnings when this file is built they this is a common origin for problems).
.. Firstly, let's add technical info as a sidebar and allow text below to wrap around it. This list is a work in
progress, please help us improve it. We use *definition lists* of ReST_ to make this readable.
.. sidebar:: Software Technical Information
Name
Thermodynamic Force Calculator for Abrupt AdResS
Language
bash
Python 2.7
Licence
MD Simulation:
See GROMACS web page: `<http://www.gromacs.org/>`_
Analysis tools:
see VOTCA web page: `<http://www.votca.org/home>`_
Documentation Tool
Application Documentation
See GROMACS web page: `<http://www.gromacs.org/>`_
See VOCTA web page: `<http://www.votca.org/Documentation>`_
Relevant Training Material
See GROMACS web page: `<http://www.gromacs.org/>`_
See VOCTA web page: `<http://www.votca.org/tutorials>`_
.. In the next line you have the name of how this module will be referenced in the main documentation (which you can
reference, in this case, as ":ref:`example`"). You *MUST* change the reference below from "example" to something
unique otherwise you will cause cross-referencing errors. The reference must come right before the heading for the
reference to work (so don't insert a comment between).
.. _adress_tf:
################################################
Thermodynamic Force Calculator for Abrupt AdResS
################################################
.. Let's add a local table of contents to help people navigate the page
.. contents:: :local:
We introduced with the Abrupt AdResS method a new way of coupling the different simulation regions together. That is the basis for easier implementation into other codes. The implementation of smooth coupling GC- AdResS in GROMACS has several performance problems. However, the new Abrupt AdResS presents a very straight forward way to implement a new partitioning scheme, which solves two problems which affect the performance, the neighbor list search and the generic force kernel. Furthermore, we update the implementation to address this in a way that decouples the method directly from the core of any MD code, which does not hinder the performance and makes the scheme hardware independent.
Theory, application and tests see `<https://aip.scitation.org/doi/10.1063/1.5031206>`_ or `<https://arxiv.org/abs/1806.09870>`_.
The drawback of this method is that a new (as in more direct) way to calculate the thermodynamic force is needed. While the theory is still the same, the interpolation has to be adapted.
.. Add an abstract for a *general* audience here. Write a few lines that explains the "helicopter view" of why you are
creating this module. For example, you might say that "This module is a stepping stone to incorporating XXXX effects
into YYYY process, which in turn should allow ZZZZ to be simulated. If successful, this could make it possible to
produce compound AAAA while avoiding expensive process BBBB and CCCC."
.. This is an example of what a *module* for E-CAM looks like. The original source of this page (:download:`readme.rst`) contains lots of additional comments to help you create your module (and understand ReST_ syntax) so please use this as a starting point. You are free add any level of complexity you wish (within the bounds of what ReST_ can do). More general instructions for making your contribution can be found in ":ref:`contributing`".
.. Remember that for a module to be accepted into the E-CAM repository, your source code changes in the target application must pass a number of acceptance criteria:
.. * Style *(use meaningful variable names, no global variables,...)*
.. * Source code documentation *(each function should be documented with each argument explained)*
.. * Tests *(everything you add should have either unit or regression tests)*
.. * Performance *(If what you introduce has a significant computational load you should make some performance optimisation effort using an appropriate tool. You should be able to verify that your changes have not introduced unexpected performance penalties, are threadsafe if needed,...)*
Purpose of Module
_________________
.. Keep the helper text below around in your module by just adding ".. " in front of it, which turns it into a comment
.. Give a brief overview of why the module is/was being created, explaining a little of the scientific background and how it fits into the larger picture of what you want to achieve.
.. If needed you can include latex mathematics like
.. :math:`\frac{ \sum_{t=0}^{N}f(t,k) }{N}`
.. which won't show up on GitLab/GitHub but will in final online documentation.
.. If you want to add a citation, such as [CIT2009]_. Note that citations may get rearranged, e.g., to the bottom of the "page".
.. : .. [CIT2009] A citation (as often used in journals).
The new Abrupt coupling scheme introduces a density discrepancy which is very much restricted to the interface of the atomistic region and the coarse grained region. The thermodynamic force calculator in VOTCA (implemented up to version 1.3) is designed for the more smooth coupling over a larger region in space. Thus this code cannot be used for the small area of disturbance in this new scheme.
Here we present a thermodynamic force calculator for the abrupt coupling scheme. It is a mix between bash and Python and can be applied even to the older smooth coupling scheme.
.. The interface between the regions is more fluctuating and needs a more responsive thermodynamic force but it works reasonably well.
.. The second piece of the puzzle is the spatial partitioning as we showed at the ESDW8 in Berlin and as Guzman et al. (arXiv:1711.03290v1) published recently it is possible to use a spatial partitioning for GC-AdResS.
Background Information
______________________
.. Keep the helper text below around in your module by just adding ".. " in front of it, which turns it into a comment
Abrupt AdResS presents a very straight forward way to implement a new partitioning scheme. The drawback is that the particles mix in a very narrow region in space. The distance of the molecules at that interface can be too close, which has to be compensated via a force capping. However the flux of particles at that interface is fast, which leads to a rather localized discrepancy in the density.
The thermodynamic force calculator in VOTCA (up to version 1.3) is designed for a smooth and not very much disturbed region in space. Thus a new code to calculate the thermodynamic force was needed. The thermodynamic force is calculated by calculating the gradient of the density in a specific region in space. Thus any code taking this into account can be used. For the detailed discussion of the role and the basic principles behind this force see `<https://aip.scitation.org/doi/10.1063/1.5031206>`_ or `<https://arxiv.org/abs/1806.09870>`_.
The code provided here is designed for easy adjustable and can be used on different computer architectures (knowledge of bash is of an advantage).
Building and Running
____________________
.. Keep the helper text below around in your module by just adding ".. " in front of it, which turns it into a comment
These three scripts present one way to calculate the thermodynamics force.
::
:download:`TF calculation script: <./TF_calc_water_xplit_sphere.sh>`
:download:`Density interpolation script: <./smooth_dens.sh>`
:download:`Script to call TF calculation: <./run_tf_water_xplsit_sphere.sh>`
The central script is *smooth_dens.sh*. This is a Python 2.7 script which interpolates the density and generates the gradient of the density and provides the force as an ascii table.
*TF_calc_water_xplit_sphere.sh* is controlling the MD run, and which region to interpolate, and builds the tables needed. The commands and the options used are described in `<http://www.gromacs.org/>`_ or if the spherical geometry for AdResS is used also here: `<http://www.votca.org/Documentation>`_.
*run_tf_water_xplsit_sphere.sh* provides some possible input scenarios for *TF_calc_water_xplit_sphere.sh*.
To run *run_tf_water_xplsit_sphere.sh* one has to first enter the correct region, which should be interpolated in *TF_calc_water_xplit_sphere.sh*.
::
rmin = starting point along a distance with rref as origin
rmax = end point of the interpolation region
rbox = maximal box size (xsplit: x direction; sphere the maximal radius)
rref = is the point defined in the GROMACS input file
lc = defines the binning along the x direction or the radius
prefac = is basically a weighting on the thermodynamic force (small: more iteration, but more careful approach of the target density)
Note of caution: in *run_tf_water_xplsit_sphere.sh* and *TF_calc_water_xplit_sphere.sh* the GROMACS and VOTCA version used have to be specifically sourced. Then select which option in *run_tf_water_xplsit_sphere.sh* you want to use and comment the other out and execute:
::
for a new run without a thermodynamic force to start with:
bash run_tf_water_xplsit_sphere.sh 1 20 1
for a start from an existing thermodynamic force:
bash run_tf_water_xplsit_sphere.sh 21 20 2
Source Code
___________
.. Notice the syntax of a URL reference below `Text <URL>`_
:download:`TF calculation script: <./TF_calc_water_xplit_sphere.sh>`
:download:`Density interpolation script: <./smooth_dens.sh>`
:download:`Script to call TF calculation: <./run_tf_water_xplsit_sphere.sh>`
#!/bin/bash
source /<put your directory path here>/gromacs_abrupt_adress/bin/GMXRC
export GMX_MAXCONSTRWARN=-1
i=$1
j=$2
case $3 in
[1]*)
#initial run: first density without thermodyn. Force
rm F_th_* tabletf* *.dens.s* dens_mix_* \#* *.dens_s*.xvg
./TF_calc_water_xplit_sphere.sh 0 dens.REF.out CALC.dens.out s0 1
#thermodyn. Force iterations from scratch (!!!) from $i (Start) to $j (end)
for z in `seq $i $j`; do
#sphere choice
./TF_calc_water_xplit_sphere.sh 2 dens.REF.out CALC.dens.out s$z 1
# xsplit choice
./TF_calc_water_xplit_sphere.sh 2 dens.REF.out CALC.dens.out s$z 2
done
#The scripts prints out the forces and densities for each iteration.
;;
[2]*)
#This is the restart from iteration $i to $j from an existing(!!!) thermodyn. Force.
for z in `seq $i $j`; do
#sphere choice
./TF_calc_water_xplit_sphere.sh 2 dens.REF.out CALC.dens.out s$z 1
# xsplit choice
./TF_calc_water_xplit_sphere.sh 2 dens.REF.out CALC.dens.out s$z 2
done
;;
esac
#!/usr/bin/env python
import sys
from numpy import *
from scipy.signal import savgol_filter
from scipy import integrate
from scipy import interpolate
import matplotlib.pyplot as plt
# only hybrid region and a little around
f = file('dens_mix_hy', 'r')
f_x = file('xnew_ntrpl', 'r')
tgt = loadtxt(f, usecols=(0,1,2))
tgt_x = loadtxt(f_x, usecols=(0))
r = tgt[:,0]
dr = tgt[:,1]
dtr = tgt[:,1]
xnew = tgt_x[1]
prefac = sys.argv[1]
diff_r = sys.argv[2]
T = sys.argv[3]
# kb*T in kJ/(mol K)
#kbT = kb*T*N/1000
kb = 1.38068452e-23
N = 6.022141e23
kbT = -1.0*0.00831451*T
rho_zero = tgt[:,2]
avg_rho = average(rho_zero)
#compressability: 4.5e-5
kappa = 4.5e-5
#prefactor: M(mol)/(rho**2 compress)
#M(CAT)=107.077
#M(ANN)=35.45300
# we set M=1 !!!
#==================================
#prepot = 1.0/(avg_rho*kappa)
#prefac = 1.0/(square(avg_rho)*kappa)
print "avg_rho",avg_rho
print "prefac", prefac
tck = interpolate.splrep(r, dr)
dr_s = interpolate.splev(r, tck, der=0)
dr_d = interpolate.splev(r, tck, der=1)
dr_p = dr_s*kbT*prefac
force_r = dr_d*prefac
SPL = column_stack((r, dr_s))
savetxt('dens_spline', SPL)
DAT = column_stack((r, force_r))
savetxt('dens_smooth', DAT)
POT = column_stack((r, dr_p))
savetxt('pot_smooth', POT)
plt.plot(r, dr)
plt.plot(r, dr_s)
#plt.show()
plt.plot(r,dr_p)
plt.plot(r,force_r)
#plt.show()
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment