Added user-manifold package.

This commit is contained in:
Stefan Paquay
2016-03-02 18:04:11 -05:00
parent 9b33b04183
commit 889fda57fe
80 changed files with 27941 additions and 0 deletions

61
doc/fix_manifoldforce.txt Normal file
View File

@ -0,0 +1,61 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
fix manifoldforce command :h3
[Syntax:]
fix ID group-ID manifoldforce manifold manifold-args ... :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
manifold = name of the manifold :l
manifold-args = parameters for the manifold :l
:ule
[Examples:]
fix constrain all manifoldforce sphere 5.0
[Description:]
This fix subtracts each time step from the force the component along the normal of the specified manifold.
This can be used in combination with "minimize"_minimize.html to remove overlap between particles while
keeping them (roughly) constrained to the given manifold, e.g. to set up a run with "fix nve/manifold/rattle"_fix_nve_manifold_rattle.html.
I have found that only {hftn} and {quickmin} with a very small time step perform adequately though.
:line
[Restart, fix_modify, output, run start/stop, minimize info:]
No information about this fix is written to "binary restart
files"_restart.html. None of the "fix_modify"_fix_modify.html options
are relevant to this fix. No global or per-atom quantities are stored
by this fix for access by various "output
commands"_Section_howto.html#howto_15. No parameter of this fix can
be used with the {start/stop} keywords of the "run"_run.html command.
This fix is invoked during "energy minimization"_minimize.html.
:line
[Restrictions:]
This fix is part of the USER-MANIFOLD package. It is only enabled if LAMMPS
was built with that package. See the "Making LAMMPS"_Section_start.html#start_3
section for more info.
Only use this with {min_style hftn} or {min_style quickmin}. If not, the constraints
will not be satisfied very well at all. A warning is generated if the {min_style} is
incompatible but no error.
:line
[Related commands:]
"fix nve/manifold/rattle"_fix_nve_manifold_rattle.html, "fix nvt/manifold/rattle"_fix_nvt_manifold_rattle.html

View File

@ -0,0 +1,100 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
fix nve/manifold/rattle command :h3
[Syntax:]
fix ID group-ID nve/manifold/rattle tol maxit manifold manifold-args keyword value ... :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
nve/manifold/rattle = style name of this fix command :l
tol = tolerance to which Newton iteration must converge :l
maxit = maximum number of iterations to perform :l
manifold = name of the manifold :l
manifold-args = parameters for the manifold :l
one or more keyword/value pairs may be appended :l
keyword = {every}
{every} values = N
N = print info about iteration every N steps. N = 0 means no output :pre
:ule
[Examples:]
fix 1 all nve/manifold/rattle 1e-4 10 sphere 5.0
fix step all nve/manifold/rattle 1e-8 100 ellipsoid 2.5 2.5 5.0 every 25 :pre
[Description:]
Perform constant NVE integration to update position and velocity for
atoms constrained to a curved surface (manifold) in the group each timestep. The constraint
is handled by RATTLE "(Andersen)"_#Andersen written out for the special case of
single-particle constraints as explained in "(Paquay)"_#Paquay.
V is volume; E is energy. This way, the dynamics of particles constrained to
curved surfaces can be studied. If combined with "fix langevin"_fix_langevin.html, this generates
Brownian motion of particles constrained to a curved surface. For a list of currently supported
manifolds and their parameters, see "manifolds"_user_manifolds.html.
Note that the particles must initially be close to the manifold in question. If not, RATTLE will
not be able to iterate until the constraint is satisfied, and an error is generated. For simple
manifolds this can be achieved with {region} and {create_atoms} commands, but for more complex
surfaces it might be more useful to write a script.
The manifold args may be equal-style variables, like so:
variable R equal "ramp(5.0,3.0)"
fix shrink_sphere all nve/manifold/rattle 1e-4 10 sphere v_R :pre
In this case, the manifold parameter will change in time according to the variable.
This is not a problem for the time integrator as long as the change of the manifold is slow with respect to the dynamics of the particles.
Note that if the manifold has to exert work on the particles because of these changes, the total energy might not be conserved.
:line
[Restart, fix_modify, output, run start/stop, minimize info:]
No information about this fix is written to "binary restart
files"_restart.html. None of the "fix_modify"_fix_modify.html options
are relevant to this fix. No global or per-atom quantities are stored
by this fix for access by various "output
commands"_Section_howto.html#howto_15. No parameter of this fix can
be used with the {start/stop} keywords of the "run"_run.html command.
This fix is not invoked during "energy minimization"_minimize.html.
:line
[Restrictions:]
This fix is part of the USER-MANIFOLD package. It is only enabled if LAMMPS
was built with that package. See the "Making LAMMPS"_Section_start.html#start_3
section for more info.
Only use this with {min_style hftn} or {min_style quickmin}. If not, the constraints
will not be satisfied very well at all. A warning is generated if the {min_style} is
incompatible but no error.
:line
[Related commands:]
"fix nvt/manifold/rattle"_fix_nvt_manifold_rattle.html, "fix manifoldforce"_fix_manifoldforce.html
[Default:] every = 0, tchain = 3
:line
:link(Andersen)
[(Andersen)] Andersen, J. Comp. Phys. 52, 24, (1983).
:link(Paquay)
[(Paquay)] Paquay and Kusters, Biophys. J., 110, ???, (2016), to be published,
preprint available at "arXiv:1411.3019"_http://arxiv.org/abs/1411.3019/.

View File

@ -0,0 +1,85 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
fix nvt/manifold/rattle command :h3
[Syntax:]
fix ID group-ID nvt/manifold/rattle tol maxit manifold manifold-args keyword value ... :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
nvt/manifold/rattle = style name of this fix command :l
tol = tolerance to which Newton iteration must converge :l
maxit = maximum number of iterations to perform :l
manifold = name of the manifold :l
manifold-args = parameters for the manifold :l
one or more keyword/value pairs may be appended :l
keyword = {temp} or {tchain} or {every}
{temp} values = Tstart Tstop Tdamp
Tstart, Tstop = external temperature at start/end of run
Tdamp = temperature damping parameter (time units)
{tchain} value = N
N = length of thermostat chain (1 = single thermostat)
{every} value = N
N = print info about iteration every N steps. N = 0 means no output :pre
:ule
[Examples:]
fix 1 all nvt/manifold/rattle 1e-4 10 cylinder 3.0 temp 1.0 1.0 10.0
[Description:]
This fix combines the RATTLE-based "(Andersen)"_#Andersen time integrator of "fix nve/manifold/rattle"_fix_nve_manifold_rattle.html "(Paquay)"_#Paquay with a Nose-Hoover-chain thermostat to sample the
canonical ensemble of particles constrained to a curved surface (manifold). This sampling does suffer from discretization bias of O(dt).
For a list of currently supported manifolds and their parameters, see "manifolds"_user_manifolds.html
:line
[Restart, fix_modify, output, run start/stop, minimize info:]
No information about this fix is written to "binary restart
files"_restart.html. None of the "fix_modify"_fix_modify.html options
are relevant to this fix. No global or per-atom quantities are stored
by this fix for access by various "output
commands"_Section_howto.html#howto_15. No parameter of this fix can
be used with the {start/stop} keywords of the "run"_run.html command.
This fix is not invoked during "energy minimization"_minimize.html.
:line
[Restrictions:]
This fix is part of the USER-MANIFOLD package. It is only enabled if LAMMPS
was built with that package. See the "Making LAMMPS"_Section_start.html#start_3
section for more info.
Only use this with {min_style hftn} or {min_style quickmin}. If not, the constraints
will not be satisfied very well at all. A warning is generated if the {min_style} is
incompatible but no error.
:line
[Related commands:]
"fix nve/manifold/rattle"_fix_nvt_manifold_rattle.html, "fix manifoldforce"_fix_manifoldforce.html
[Default:] every = 0
:line
:link(Andersen)
[(Andersen)] Andersen, J. Comp. Phys. 52, 24, (1983).
:link(Paquay)
[(Paquay)] Paquay and Kusters, Biophys. J., 110, ???, (2016), to be published,
preprint available at "arXiv:1411.3019"_http://arxiv.org/abs/1411.3019/.

30
doc/user_manifolds.txt Normal file
View File

@ -0,0 +1,30 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
[Description:]
Below is a list of currently supported manifolds, their parameters and a short description of them.
The parameters listed here are in the same order as they should be passed to the relevant fixes.
{manifold} || {parameters} || {equation} || {description}
cylinder || R || x^2 + y^2 - R^2 = 0 || Cylinder along z-axis, axis going through (0,0,0)
cylinder_dent || R l a || x^2 + y^2 - r(z)^2 = 0, r(x) = R if |z| > l, r(z) = R - a*(1 + cos(z/l))/2 otherwise || A cylinder with a dent around z = 0
dumbbell || a A B c || -( x^2 + y^2 ) * (a^2 - z^2/c^2) * ( 1 + (A*sin(B*z^2))^4) = 0 || A dumbbell ||
ellipsoid || a b c || (x/a)^2 + (y/b)^2 + (z/c)^2 = 0 || An ellipsoid
plane || a b c x0 y0 z0 || a*(x-x0) + b*(y-y0) + c*(z-z0) = 0 || A plane with normal (a,b,c) going through point (x0,y0,z0)
plane_wiggle || a w || z - a*sin(w*x) = 0 || A plane with a sinusoidal modulation on z along x.
sphere || R || x^2 + y^2 + z^2 - R^2 = 0 || A sphere of radius R
supersphere || R q || |x|^q + |y|^q + |z|^q - R^q = 0 || A supersphere of hyperradius R
spine || a, A, B, B2, c || -(x^2 + y^2)*(a^2 - z^2/f(z)^2)*(1 + (A*sin(g(z)*z^2))^4), f(z) = c if z > 0, 1 otherwise; g(z) = B if z > 0, B2 otherwise || An approximation to a dendtritic spine
thylakoid || wB LB lB || Various, see "(Paquay)"_#Paquay || A model grana thylakoid consisting of two block-like compartments connected by a bridge of width wB, length LB and taper length lB
torus || R r || (R - sqrt( x^2 + y^2 ) )^2 + z^2 - r^2 || A torus with large radius R and small radius r, centered on (0,0,0) :tb(s=||)
:link(Paquay)
[(Paquay)] Paquay and Kusters, Biophys. J., 110, ???, (2016), to be published,
preprint available at "arXiv:1411.3019"_http://arxiv.org/abs/1411.3019/.

View File

@ -0,0 +1,4 @@
[Dolphin]
Timestamp=2016,3,2,16,8,35
Version=3
ViewMode=1

View File

@ -0,0 +1 @@
This example uses fix nve/manifold/rattle and a Langevin thermostat to generate Brownian motion on a few manifolds. A plotting script is included to plot the mean squared displacement against theoretical expectation.

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,30 @@
dimension 3
units lj
boundary f f p
read_data msd.cyl.data
mass 1 1.0
pair_style none
atom_modify sort 0 1.0
variable R equal 5
fix step all nve/manifold/rattle 1e-10 100 cylinder $R
fix temp all langevin 1.0 1.0 1.0 12321 gjf yes
velocity all create 1.0 1283
run 25000
compute dx2 all msd
variable D equal "1.0"
variable t equal time
variable my_msd equal "2*v_D*(v_t-125)+ 2*v_R*v_R*(1.0 - exp(-(v_t-125)*v_D/(v_R*v_R)) )"
variable msd equal c_dx2[4]
dump traj all custom 100 msd.cyl.gz id type x y zu
fix msd_out all print 250 "$t ${msd} ${my_msd}" file msd.cyl.dat screen no
thermo_style custom time step pe ke etotal temp c_dx2[4] v_my_msd
thermo 1000
run 25000

Binary file not shown.

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,29 @@
dimension 3
units lj
boundary p p f
read_data msd.plane.data
mass 1 1.0
pair_style none
atom_modify sort 0 1.0
variable R equal 8
fix step all nve/manifold/rattle 1e-4 4 plane 0 0 1 0 0 0
fix temp all langevin 1.0 1.0 1.0 12321 gjf yes
velocity all create 1.0 1283
run 25000
compute dx2 all msd
variable D equal "1.0"
variable t equal time
variable my_msd equal "4*(v_t - 125)*v_D"
variable msd equal c_dx2[4]
fix msd_out all print 250 "$t ${msd} ${my_msd}" file msd.plane.dat screen no
thermo_style custom time step pe ke etotal temp c_dx2[4] v_my_msd
thermo 1000
run 25000

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,27 @@
dimension 3
units lj
boundary f f f
read_data msd.sphere.data
mass 1 1.0
pair_style none
atom_modify sort 0 1.0
variable R equal 8
fix step all nve/manifold/rattle 1e-4 4 sphere $R
fix temp all langevin 1.0 1.0 1.0 12321 gjf yes
velocity all create 1.0 1283
run 25000
compute dx2 all msd
variable D equal "1.0"
variable t equal time
variable my_msd equal "2*v_R*v_R*(1.0 - exp(-2*(v_t-125)*v_D/(v_R*v_R)) )"
variable msd equal c_dx2[4]
fix msd_out all print 250 "$t ${msd} ${my_msd}" file msd.sphere.dat screen no
thermo_style custom time step pe ke etotal temp c_dx2[4] v_my_msd
thermo 1000
run 25000

Binary file not shown.

View File

@ -0,0 +1,63 @@
# Plots the observed and expected MSDs.
#
set term epslatex standalone size 12cm,8cm color dashed
set output 'msd_plot2.tex'
msd_s(t,R,D) = 2*R*R*( 1 - exp(-2*D*t/(R*R)) )
msd_c(t,R,D) = 2*D*t + 2*R*R*( 1 - exp(-D*t/(R*R)) )
msd_p(t,D) = 4*D*t
set key bottom right box opaque maxrow 3 samplen 1 width -5
set title 'Mean Squared Displacements on curved surfaces'
set xlabel 'Time / $\tau$' offset 0,0.5
set ylabel 'Mean squared displacement $\left<\delta x^2\right>$'
set multiplot
set size 1,1
set origin 0,0
set xrange[0:130]
set yrange[-80:520]
set style line 1 lc rgb '#000000' lw 6.0 dt 2
set style line 2 lc rgb '#000000' lw 6.0 dt 3
set style line 3 lc rgb '#000000' lw 3.0 dt 1
set grid
plot 'msd.plane.dat' u ($1-125):2 w p pt 4 ps 2.0 lw 2 lc rgb '#00AA00' ti 'Plane, data' \
, 'msd.cyl.dat' u ($1-125):2 w p pt 6 ps 2.0 lw 2 lc rgb '#AA0000' ti 'Cylinder, data' \
, 'msd.sphere.dat' u ($1-125):2 w p pt 8 ps 2.0 lw 2 lc rgb '#0000CC' ti 'Sphere, data' \
, msd_p(x, 1.0) w l ls 1 ti 'Plane, theory' \
, msd_c(x, 5, 1.0) w l ls 2 ti 'Cylinder, theory' \
, msd_s(x, 8, 1.0) w l ls 3 ti 'Sphere, theory'
set origin 0.12,0.46
set size 0.4,0.42
# Hack to remove old grid from inset:
set object 1 rectangle from graph 0,0 to graph 1,1 fillcolor rgb "white" behind
unset key
set grid front
set title 'Short time behaviour' offset 0,-0.8
set ylabel ''
set xrange[0:10]
set yrange[0:40]
set ytics 10
set xtics 2
plot 'msd.plane.dat' u ($1-125):2 w p pt 4 ps 2.0 lw 2 lc rgb '#00AA00' ti 'Plane, data' \
, 'msd.cyl.dat' u ($1-125):2 w p pt 6 ps 2.0 lw 2 lc rgb '#AA0000' ti 'Cylinder, data' \
, 'msd.sphere.dat' u ($1-125):2 w p pt 8 ps 2.0 lw 2 lc rgb '#0000CC' ti 'Sphere, data' \
, msd_p(x, 1.0) w l ls 1 ti 'plane, theory' \
, msd_c(x, 5, 1.0) w l ls 2 ti 'cylinder, theory' \
, msd_s(x, 8, 1.0) w l ls 3 ti 'sphere, theory'
unset multiplot

View File

@ -0,0 +1,4 @@
[Dolphin]
Timestamp=2016,3,2,16,6,15
Version=3
ViewMode=1

View File

@ -0,0 +1,3 @@
This example runs an NVE simulation on various manifolds and checks for energy conservation.
The script "energy.sh" automates the running of LAMMPS for the different manifolds, and if gnuplot is installed, it makes a nice plot.

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,3 @@
boundary f f p
read_data cylinder.data
variable manifold_params string "6.0"

View File

@ -0,0 +1,28 @@
units lj
dimension 3
include ${manifold}.setup
fix step all nve/manifold/rattle 1e-6 10 ${manifold} ${manifold_params}
pair_style lj/smooth/linear 2.5
pair_coeff * * 1.0 1.0
pair_modify shift yes
# Rebalance:
fix load all balance 1000 1.0 shift "xyz" 10 1.05
timestep 0.0005
variable t equal time
variable U equal pe
variable K equal ke
variable E equal "v_U + v_K"
fix therm all print 25000 "$t $U $K $E" file thermo.${manifold}.dat screen no
dump traj all custom 1000000 dump.${manifold}.gz id type x y z
thermo_style custom time step temp pe ke etotal
thermo 10000
run ${STEPS}

View File

@ -0,0 +1,27 @@
units lj
dimension 2
read_data plane.data
fix step all nve
fix twod all enforce2d
pair_style lj/smooth/linear 2.5
pair_coeff * * 1.0 1.0
pair_modify shift yes
timestep 0.0005
variable t equal time
variable U equal pe
variable K equal ke
variable E equal "v_U + v_K"
fix therm all print 25000 "$t $U $K $E" file thermo.plane.dat screen no
dump traj all custom 1000000 dump.plane.gz id type x y z
thermo_style custom time step temp pe ke etotal
thermo 10000
run ${STEPS}

View File

@ -0,0 +1,39 @@
#!/bin/bash
#
# Runs the various tests.
#
# Set this to your lammps executable
LMP="mpirun -np 2 lmp"
NO_GNUPLOT=$( command -v gnuplot >/dev/null 2>&1 )
NSTEPS=100000000
if [ $# -ge 1 ]; then
NSTEPS=$1
fi
for MANIFOLD in torus sphere cylinder
do
echo "+=== Running example on "$MANIFOLD" ===+"
DATE=$(date | awk '{ print $3$2substr($6,3) }')
LOG=$MANIFOLD".energy."$DATE".log"
$LMP -in energy.in -var manifold $MANIFOLD -log $MANIFOLD".energy.log" -var STEPS $NSTEPS
echo "+======================================+"
echo ""
done
# 2D plane for verification:
$LMP -in energy.plane.in -log energy.plane.log -var STEPS $NSTEPS
if [ $NO_GNUPLOT ]
then
echo "No Gnuplot found, not plotting."
exit 0
fi
E0S=$( head -n 2 thermo.sphere.dat | tail -n 1 | awk '{ print $4 }' )
E0C=$( head -n 2 thermo.cylinder.dat | tail -n 1 | awk '{ print $4 }' )
E0T=$( head -n 2 thermo.torus.dat | tail -n 1 | awk '{ print $4 }' )
E0P=$( head -n 2 thermo.plane.dat | tail -n 1 | awk '{ print $4 }' )
echo "Plotting using Gnuplot"
gnuplot -e "E0S="$E0S -e "E0C="$E0C -e "E0T="$E0T -e "E0P="$E0P plot_energies.gpl

Binary file not shown.

Binary file not shown.

After

Width:  |  Height:  |  Size: 84 KiB

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,78 @@
# Gnuplot script to plot the energies.
#
#set term pngcairo size 800,600 solid color
#set output 'energy_conservation.png'
set term epslatex size 14cm,10cm standalone color
set output 'energy_conservation.tex'
E0S = 0.263720171785176
E0C = 0.680969247210183
E0T = 0.393139981206627
E0P = -1.16336568585219
E1S = 0.168915271276422
E1C = 0.615216219988708
E1T = 0.318731094092091
E1P = -1.29808304730092
set grid
set xlabel 'Time (LJ units)'
set ylabel 'Relative energy deviation $(\mathcal{H}(t) - \mathcal{H}(0))/\mathcal{H}(0)$'
set title 'Energy drift for RATTLE-M on various surfaces'
set key opaque box top left
set log x
col_S = '#CC0000'
col_C = '#00CC00'
col_T = '#0000CC'
col_P = '#000000'
set multiplot
set xrange[10:5e6]
set yrange[-0.02:0.02]
set ytics ( "$-2.0~10^{-2}$" -2e-2, "$-1.5~10^{-2}$" -1.5e-2, "$-1.0~10^{-2}$" -1e-2 \
, "$-0.5~10^{-2}$" -0.5e-2, "$0$" 0, "$0.5~10^{-2}$" 0.5e-2, "$1.0~10^{-2}$" 1e-2 \
, "$1.5~10^{-2}$" 1.5e-2, "$2.0~10^{-2}$" 2e-2 )
plot 'long_run_lj_cut/thermo.sphere.dat' u 1:($4/E1S-1) w l lw 2 lc rgb col_S ti 'S, Truncated' \
, 'long_run_lj_cut/thermo.cylinder.dat' u 1:($4/E1C-1) w l lw 2 lc rgb col_C ti 'C, Truncated' \
, 'long_run_lj_cut/thermo.torus.dat' u 1:($4/E1T-1) w l lw 2 lc rgb col_T ti 'T, Truncated' \
, 'long_run_lj_cut/thermo.plane.dat' u 1:($4/E1P-1) w l lw 2 lc rgb col_P ti 'P, Truncated' \
, 'thermo.sphere.dat' u 1:($4/E0S-1) w p pt 2 ps 1 lc rgb col_S ti 'S, Smoothed' \
, 'thermo.cylinder.dat' u 1:($4/E0C-1) w p pt 4 ps 1 lc rgb col_C ti 'C, Smoothed' \
, 'thermo.torus.dat' u 1:($4/E0T-1) w p pt 6 ps 1 lc rgb col_T ti 'T, Smoothed' \
, 'thermo.plane.dat' u 1:($4/E0P-1) w p pt 8 ps 1 lc rgb col_P ti 'P, Smoothed'
set origin 0.285,0.106
set size 0.55,0.4
set object 1 rectangle from graph 0,0 to graph 1,1 fillcolor rgb "white" behind
set xlabel ''
set ylabel ''
set title ''
set xrange[10:500000]
set yrange[-1.5e-5:1.5e-5]
set ytics
set format y "%0.1t 10^{%T}"
set grid
unset key
#set ytics ( "$-1.5~10^{-5}$" -1.5e-5, "$-1.0~10^{-5}$" -1e-5, "$-0.5~10^{-5}$" -0.5e-5, \
# "$0$" 0, "$0.5~10^{-5}$" 0.5e-5, "$1.0~10^{-5}$" 1.0e-5, "$1.5~10^{-5}$" 1.5e-5 )
set ytics ( "$-1.8~10^{-5}$" -1.8e-5, "$-1.2~10^{-5}$" -1.2e-5, "$-0.6~10^{-5}$" -0.6e-5, \
"$0$" 0, "$0.6~10^{-5}$" 0.6e-5, "$1.2~10^{-5}$" 1.2e-5, "$1.8~10^{-5}$" 1.8e-5 )
set title 'Smoothed only' offset 0,-0.5
plot 'thermo.sphere.dat' u 1:($4/E0S - 1) w p pt 2 ps 1 lc rgb col_S \
, 'thermo.cylinder.dat' u 1:($4/E0C - 1) w p pt 4 ps 1 lc rgb col_C \
, 'thermo.torus.dat' u 1:($4/E0T - 1) w p pt 6 ps 1 lc rgb col_T \
, 'thermo.plane.dat' u 1:($4/E0P - 1) w p pt 8 ps 1 lc rgb col_P
unset multiplot

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,3 @@
boundary f f f
read_data sphere.data
variable manifold_params string "10.0"

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,3 @@
boundary f f f
read_data torus.data
variable manifold_params string "10.0 4.0"

View File

@ -0,0 +1,4 @@
[Dolphin]
Timestamp=2016,3,2,16,13,51
Version=3
ViewMode=1

View File

@ -0,0 +1 @@
This example simulates particles on a complex manifold that consists of pieces of cylinder, plane and sphere. By keeping track of when the blue particles reach the other compartment, information about mean escape times can be obtained.

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,33 @@
units lj
dimension 3
boundary f f p
pair_style lj/cut 1.1225
read_data thylakoid.data
pair_coeff * * 1.0 1.0
pair_modify shift yes
# Makes about 100 particles type 2:
region make_two sphere -29 0 0 6.25
region everywhere block -25 25 -25 25 -25 25
group two region make_two
set group two type 2
# And deletes some of type 1 to allow more movement:
group one type 1
set group one type 1
delete_atoms porosity everywhere 0.2 123
# Start run
variable wB equal 5.0
variable LB equal 10.0
variable lB equal 3.0
fix step all nve/manifold/rattle 1e-4 10 thylakoid ${wB} ${LB} ${lB}
thermo_style custom time step pe ke etotal temp
thermo 2500
dump traj all custom 10000 thylakoid.dump id type x y z
run 25000000

View File

@ -0,0 +1,4 @@
[Dolphin]
Timestamp=2016,3,2,16,11,43
Version=3
ViewMode=1

View File

@ -0,0 +1 @@
This example simulates some coarse-grained capsomere model attached to a slowly shrinking spherical surface.

View File

@ -0,0 +1,704 @@
LAMMPS data file with capsid model
288 atoms
4 atom types
216 bonds
3 bond types
144 angles
1 angle types
-8 8 xlo xhi
-8 8 ylo yhi
-8 8 zlo zhi
Masses
1 1.0
2 1.1.728
3 2.744
4 1.0
Pair Coeffs # lj/cut
1 1 1 1.1225
2 1 1.35 4
3 1 2.1 4
4 1 1 1.1225
PairIJ Coeffs # lj/cut
1 1 1 1 1.1225
1 2 1 1.175 1.31894
1 3 1 1.55 1.73988
1 4 1 1 1.1225
2 2 1 1.35 4
2 3 1 1.725 1.93631
2 4 1 1.175 1.31894
3 3 1 2.1 4
3 4 1 1.55 1.73988
4 4 1 1 1.1225
Bond Coeffs # harmonic
1 50 0.659469
2 50 0.855906
3 50 0.659469
Atoms
1 1 1 -0.701084 3.93343 0.19124
2 1 2 -0.816669 4.58193 0.222769
3 1 3 -0.966685 5.42359 0.263689
4 1 4 -1.08227 6.07209 0.295219
5 2 1 0.0618129 0.850751 3.90799
6 2 2 0.0720038 0.991012 4.55229
7 2 3 0.0852304 1.17305 5.38851
8 2 4 0.0954213 1.31331 6.03281
9 3 1 3.30064 -0.413296 -2.22148
10 3 2 3.8448 -0.481435 -2.58773
11 3 3 4.55106 -0.569871 -3.06307
12 3 4 5.09523 -0.638009 -3.42932
13 4 1 -1.59122 -1.3139 -3.42661
14 4 2 -1.85356 -1.53052 -3.99155
15 4 3 -2.19405 -1.81166 -4.72477
16 4 4 -2.45639 -2.02828 -5.2897
17 5 1 -3.9047 -0.80606 -0.32184
18 5 2 -4.54846 -0.938953 -0.374901
19 5 3 -5.38397 -1.11143 -0.443767
20 5 4 -6.02773 -1.24432 -0.496828
21 6 1 1.01166 3.7639 -0.899789
22 6 2 1.17845 4.38444 -1.04813
23 6 3 1.39492 5.18983 -1.24067
24 6 4 1.56171 5.81037 -1.38901
25 7 1 -1.05538 3.74967 -0.908937
26 7 2 -1.22938 4.36786 -1.05879
27 7 3 -1.45521 5.1702 -1.25328
28 7 4 -1.62921 5.7884 -1.40314
29 8 1 -1.51588 -2.46782 2.75898
30 8 2 -1.76579 -2.87469 3.21385
31 8 3 -2.09016 -3.40274 3.8042
32 8 4 -2.34007 -3.8096 4.25907
33 9 1 -2.25141 0.367776 3.28571
34 9 2 -2.62259 0.42841 3.82742
35 9 3 -3.10434 0.507105 4.53048
36 9 4 -3.47552 0.56774 5.07219
37 10 1 1.67684 1.52227 3.2971
38 10 2 1.9533 1.77324 3.84069
39 10 3 2.31211 2.09897 4.54619
40 10 4 2.58856 2.34994 5.08978
41 11 1 -3.25052 2.13047 -0.94615
42 11 2 -3.78643 2.48171 -1.10214
43 11 3 -4.48196 2.93758 -1.30459
44 11 4 -5.01787 3.28883 -1.46058
45 12 1 0.0676891 -0.567005 3.95903
46 12 2 0.0788488 -0.660486 4.61174
47 12 3 0.0933327 -0.781812 5.45888
48 12 4 0.104492 -0.875292 6.1116
49 13 1 1.59585 -1.27299 3.43988
50 13 2 1.85895 -1.48286 4.007
51 13 3 2.20043 -1.75525 4.74306
52 13 4 2.46353 -1.96512 5.31018
53 14 1 -1.40854 2.22169 -3.01332
54 14 2 -1.64076 2.58797 -3.51012
55 14 3 -1.94215 3.06336 -4.1549
56 14 4 -2.17438 3.42964 -4.6517
57 15 1 1.7153 -2.87555 -2.18838
58 15 2 1.99809 -3.34963 -2.54917
59 15 3 2.36513 -3.96493 -3.01743
60 15 4 2.64792 -4.43901 -3.37822
61 16 1 -2.22107 -2.07058 -2.60376
62 16 2 -2.58725 -2.41195 -3.03303
63 16 3 -3.06251 -2.85501 -3.59017
64 16 4 -3.42869 -3.19638 -4.01945
65 17 1 2.2192 1.29232 -3.06677
66 17 2 2.58507 1.50538 -3.57238
67 17 3 3.05993 1.7819 -4.2286
68 17 4 3.4258 1.99496 -4.73421
69 18 1 -2.03962 2.69954 -2.13364
70 18 2 -2.37588 3.14461 -2.48541
71 18 3 -2.81231 3.72225 -2.94196
72 18 4 -3.14858 4.16731 -3.29373
73 19 1 -3.49737 -0.254473 1.92449
74 19 2 -4.07397 -0.296428 2.24178
75 19 3 -4.82232 -0.350879 2.65357
76 19 4 -5.39892 -0.392833 2.97086
77 20 1 3.50808 -1.84615 0.533931
78 20 2 4.08645 -2.15052 0.621959
79 20 3 4.8371 -2.54555 0.736207
80 20 4 5.41547 -2.84992 0.824235
81 21 1 -0.325119 3.30809 2.22505
82 21 2 -0.378721 3.85349 2.59188
83 21 3 -0.448289 4.56134 3.06799
84 21 4 -0.50189 5.10674 3.43483
85 22 1 0.662812 0.0220505 -3.94464
86 22 2 0.772088 0.0256859 -4.59498
87 22 3 0.913914 0.0304042 -5.43904
88 22 4 1.02319 0.0340396 -6.08939
89 23 1 3.23422 -1.85295 -1.45135
90 23 2 3.76744 -2.15844 -1.69062
91 23 3 4.45949 -2.55492 -2.00118
92 23 4 4.9927 -2.86041 -2.24046
93 24 1 -3.96042 -0.0236549 0.560805
94 24 2 -4.61337 -0.0275549 0.653263
95 24 3 -5.4608 -0.0326165 0.773262
96 24 4 -6.11375 -0.0365164 0.865721
97 25 1 -0.68153 -3.49389 -1.82434
98 25 2 -0.793892 -4.06992 -2.12512
99 25 3 -0.939723 -4.81753 -2.51548
100 25 4 -1.05208 -5.39356 -2.81626
101 26 1 0.487017 3.94319 0.462684
102 26 2 0.56731 4.59329 0.538966
103 26 3 0.671521 5.43704 0.63797
104 26 4 0.751814 6.08714 0.714251
105 27 1 3.01581 -0.475516 2.58434
106 27 2 3.51301 -0.553913 3.01041
107 27 3 4.15833 -0.655663 3.5634
108 27 4 4.65553 -0.73406 3.98947
109 28 1 -0.782228 -0.526685 -3.88725
110 28 2 -0.911192 -0.613518 -4.52813
111 28 3 -1.07857 -0.726216 -5.35991
112 28 4 -1.20753 -0.813049 -6.00079
113 29 1 -2.59565 -0.97452 2.88321
114 29 2 -3.02359 -1.13519 3.35856
115 29 3 -3.579 -1.34371 3.97549
116 29 4 -4.00694 -1.50438 4.45084
117 30 1 1.7976 -0.409766 -3.54975
118 30 2 2.09396 -0.477322 -4.13499
119 30 3 2.4786 -0.565003 -4.89455
120 30 4 2.77497 -0.63256 -5.47979
121 31 1 1.32145 -3.66441 -0.908779
122 31 2 1.53931 -4.26855 -1.05861
123 31 3 1.82207 -5.05265 -1.25306
124 31 4 2.03993 -5.65679 -1.40289
125 32 1 0.541532 2.8122 2.79254
126 32 2 0.630813 3.27584 3.25294
127 32 3 0.746688 3.87758 3.85048
128 32 4 0.835969 4.34122 4.31088
129 33 1 0.85239 -2.98497 2.52257
130 33 2 0.992921 -3.4771 2.93846
131 33 3 1.17531 -4.11581 3.47823
132 33 4 1.31584 -4.60794 3.89412
133 34 1 2.3452 3.20727 0.46198
134 34 2 2.73185 3.73605 0.538145
135 34 3 3.23366 4.42233 0.636998
136 34 4 3.62031 4.9511 0.713164
137 35 1 -0.487077 -3.85764 0.938818
138 35 2 -0.56738 -4.49364 1.0936
139 35 3 -0.671604 -5.31908 1.29448
140 35 4 -0.751907 -5.95508 1.44926
141 36 1 2.5761 1.56047 2.63223
142 36 2 3.00081 1.81774 3.0662
143 36 3 3.55204 2.15164 3.62944
144 36 4 3.97675 2.40892 4.06341
145 37 1 1.1532 -1.61546 -3.47281
146 37 2 1.34333 -1.8818 -4.04536
147 37 3 1.59008 -2.22747 -4.78846
148 37 4 1.78021 -2.49381 -5.36101
149 38 1 3.80441 1.23346 -0.0709251
150 38 2 4.43164 1.43681 -0.0826184
151 38 3 5.24569 1.70074 -0.0977947
152 38 4 5.87291 1.9041 -0.109488
153 39 1 -0.831367 -1.75837 3.49528
154 39 2 -0.968432 -2.04827 4.07153
155 39 3 -1.14632 -2.42452 4.81944
156 39 4 -1.28339 -2.71442 5.39569
157 40 1 -3.35775 -2.00673 -0.835774
158 40 2 -3.91134 -2.33757 -0.973565
159 40 3 -4.62982 -2.76697 -1.1524
160 40 4 -5.1834 -3.09781 -1.29019
161 41 1 -2.20432 1.30999 -3.07
162 41 2 -2.56774 1.52596 -3.57614
163 41 3 -3.03942 1.80627 -4.23305
164 41 4 -3.40284 2.02224 -4.73919
165 42 1 1.09288 -3.66612 1.16842
166 42 2 1.27306 -4.27054 1.36105
167 42 3 1.50691 -5.055 1.61106
168 42 4 1.68709 -5.65943 1.8037
169 43 1 -3.07386 0.781288 2.43741
170 43 2 -3.58064 0.910097 2.83926
171 43 3 -4.23837 1.07727 3.36081
172 43 4 -4.74515 1.20608 3.76266
173 44 1 -3.15076 -2.10328 1.28411
174 44 2 -3.67021 -2.45004 1.49582
175 44 3 -4.3444 -2.90009 1.77059
176 44 4 -4.86386 -3.24686 1.9823
177 45 1 1.83257 3.15105 -1.64698
178 45 2 2.1347 3.67056 -1.91852
179 45 3 2.52683 4.34481 -2.27093
180 45 4 2.82896 4.86431 -2.54247
181 46 1 3.03229 2.44072 0.920939
182 46 2 3.53221 2.84312 1.07277
183 46 3 4.18105 3.36537 1.26983
184 46 4 4.68097 3.76777 1.42166
185 47 1 -0.685148 -2.46174 -3.0774
186 47 2 -0.798106 -2.8676 -3.58476
187 47 3 -0.944712 -3.39436 -4.24326
188 47 4 -1.05767 -3.80022 -4.75062
189 48 1 1.86019 2.48267 2.52509
190 48 2 2.16687 2.89198 2.94139
191 48 3 2.5649 3.42321 3.4817
192 48 4 2.87159 3.83252 3.898
193 49 1 3.92074 -0.28145 -0.740678
194 49 2 4.56714 -0.327852 -0.862792
195 49 3 5.40608 -0.388076 -1.02128
196 49 4 6.05249 -0.434478 -1.14339
197 50 1 -2.4331 -3.17132 -0.150833
198 50 2 -2.83424 -3.69417 -0.1757
199 50 3 -3.35486 -4.37276 -0.207975
200 50 4 -3.756 -4.8956 -0.232842
201 51 1 -2.9946 1.69085 -2.04289
202 51 2 -3.48831 1.96961 -2.3797
203 51 3 -4.12908 2.33142 -2.81683
204 51 4 -4.62279 2.61018 -3.15364
205 52 1 0.954769 3.56359 1.54571
206 52 2 1.11218 4.15111 1.80055
207 52 3 1.31648 4.91364 2.13129
208 52 4 1.47389 5.50116 2.38613
209 53 1 -0.000383688 3.35883 -2.17216
210 53 2 -0.000446945 3.91259 -2.53028
211 53 3 -0.000529046 4.6313 -2.99507
212 53 4 -0.000592303 5.18506 -3.35318
213 54 1 -1.13409 2.01537 3.26376
214 54 2 -1.32106 2.34763 3.80185
215 54 3 -1.56373 2.77887 4.50022
216 54 4 -1.7507 3.11114 5.03831
217 55 1 -0.18323 1.16707 -3.82157
218 55 2 -0.213438 1.35949 -4.45162
219 55 3 -0.252645 1.60921 -5.26934
220 55 4 -0.282854 1.80162 -5.89939
221 56 1 -3.16152 -1.58871 -1.86569
222 56 2 -3.68275 -1.85063 -2.17328
223 56 3 -4.35924 -2.19058 -2.57249
224 56 4 -4.88048 -2.4525 -2.88008
225 57 1 0.32697 -3.16727 -2.42105
226 57 2 0.380876 -3.68945 -2.82021
227 57 3 0.45084 -4.36717 -3.33825
228 57 4 0.504746 -4.88935 -3.73741
229 58 1 -2.57648 2.36979 1.93542
230 58 2 -3.00126 2.76049 2.25451
231 58 3 -3.55256 3.26756 2.66865
232 58 4 -3.97734 3.65826 2.98773
233 59 1 -1.84228 -3.2137 1.50934
234 59 2 -2.14601 -3.74354 1.75819
235 59 3 -2.54021 -4.4312 2.08115
236 59 4 -2.84394 -4.96103 2.32999
237 60 1 2.06726 -2.51456 2.32453
238 60 2 2.40808 -2.92913 2.70777
239 60 3 2.85043 -3.46718 3.20516
240 60 4 3.19125 -3.88175 3.5884
241 61 1 3.30924 -1.38831 1.76679
242 61 2 3.85482 -1.6172 2.05807
243 61 3 4.56292 -1.91427 2.43612
244 61 4 5.10851 -2.14315 2.72741
245 62 1 0.0662913 -3.53779 1.86538
246 62 2 0.0772206 -4.12106 2.17292
247 62 3 0.0914054 -4.87806 2.57206
248 62 4 0.102335 -5.46133 2.8796
249 63 1 -0.984186 3.66087 1.27648
250 63 2 -1.14645 4.26443 1.48693
251 63 3 -1.35704 5.04777 1.76007
252 63 4 -1.5193 5.65133 1.97052
253 64 1 -2.57041 2.98875 -0.678513
254 64 2 -2.99419 3.48149 -0.790377
255 64 3 -3.54419 4.12102 -0.935563
256 64 4 -3.96797 4.61376 -1.04743
257 65 1 -0.127433 -3.91751 -0.798032
258 65 2 -0.148443 -4.56338 -0.929602
259 65 3 -0.175711 -5.40164 -1.10036
260 65 4 -0.19672 -6.04751 -1.23193
261 66 1 -3.68506 0.947699 1.23377
262 66 2 -4.29261 1.10394 1.43718
263 66 3 -5.08112 1.30673 1.70118
264 66 4 -5.68867 1.46297 1.90459
265 67 1 -2.16602 -2.87045 -1.75181
266 67 2 -2.52313 -3.3437 -2.04062
267 67 3 -2.98661 -3.95791 -2.41547
268 67 4 -3.34371 -4.43115 -2.70429
269 68 1 -0.969111 0.140347 3.87829
270 68 2 -1.12889 0.163486 4.51769
271 68 3 -1.33625 0.193517 5.34755
272 68 4 -1.49603 0.216656 5.98696
273 69 1 3.55906 0.798001 -1.64203
274 69 2 4.14584 0.929565 -1.91275
275 69 3 4.90739 1.10032 -2.2641
276 69 4 5.49417 1.23188 -2.53482
277 70 1 -3.40994 -0.202713 -2.08115
278 70 2 -3.97213 -0.236134 -2.42427
279 70 3 -4.70178 -0.279509 -2.86958
280 70 4 -5.26397 -0.31293 -3.2127
281 71 1 2.62812 -1.82176 -2.40295
282 71 2 3.06141 -2.12211 -2.79912
283 71 3 3.62377 -2.51192 -3.31329
284 71 4 4.05706 -2.81227 -3.70946
285 72 1 -3.82746 0.847616 -0.795064
286 72 2 -4.45848 0.98736 -0.926144
287 72 3 -5.27746 1.16873 -1.09627
288 72 4 -5.90849 1.30847 -1.22735
Bonds
1 1 1 2
2 2 2 3
3 3 3 4
4 1 5 6
5 2 6 7
6 3 7 8
7 1 9 10
8 2 10 11
9 3 11 12
10 1 13 14
11 2 14 15
12 3 15 16
13 1 17 18
14 2 18 19
15 3 19 20
16 1 21 22
17 2 22 23
18 3 23 24
19 1 25 26
20 2 26 27
21 3 27 28
22 1 29 30
23 2 30 31
24 3 31 32
25 1 33 34
26 2 34 35
27 3 35 36
28 1 37 38
29 2 38 39
30 3 39 40
31 1 41 42
32 2 42 43
33 3 43 44
34 1 45 46
35 2 46 47
36 3 47 48
37 1 49 50
38 2 50 51
39 3 51 52
40 1 53 54
41 2 54 55
42 3 55 56
43 1 57 58
44 2 58 59
45 3 59 60
46 1 61 62
47 2 62 63
48 3 63 64
49 1 65 66
50 2 66 67
51 3 67 68
52 1 69 70
53 2 70 71
54 3 71 72
55 1 73 74
56 2 74 75
57 3 75 76
58 1 77 78
59 2 78 79
60 3 79 80
61 1 81 82
62 2 82 83
63 3 83 84
64 1 85 86
65 2 86 87
66 3 87 88
67 1 89 90
68 2 90 91
69 3 91 92
70 1 93 94
71 2 94 95
72 3 95 96
73 1 97 98
74 2 98 99
75 3 99 100
76 1 101 102
77 2 102 103
78 3 103 104
79 1 105 106
80 2 106 107
81 3 107 108
82 1 109 110
83 2 110 111
84 3 111 112
85 1 113 114
86 2 114 115
87 3 115 116
88 1 117 118
89 2 118 119
90 3 119 120
91 1 121 122
92 2 122 123
93 3 123 124
94 1 125 126
95 2 126 127
96 3 127 128
97 1 129 130
98 2 130 131
99 3 131 132
100 1 133 134
101 2 134 135
102 3 135 136
103 1 137 138
104 2 138 139
105 3 139 140
106 1 141 142
107 2 142 143
108 3 143 144
109 1 145 146
110 2 146 147
111 3 147 148
112 1 149 150
113 2 150 151
114 3 151 152
115 1 153 154
116 2 154 155
117 3 155 156
118 1 157 158
119 2 158 159
120 3 159 160
121 1 161 162
122 2 162 163
123 3 163 164
124 1 165 166
125 2 166 167
126 3 167 168
127 1 169 170
128 2 170 171
129 3 171 172
130 1 173 174
131 2 174 175
132 3 175 176
133 1 177 178
134 2 178 179
135 3 179 180
136 1 181 182
137 2 182 183
138 3 183 184
139 1 185 186
140 2 186 187
141 3 187 188
142 1 189 190
143 2 190 191
144 3 191 192
145 1 193 194
146 2 194 195
147 3 195 196
148 1 197 198
149 2 198 199
150 3 199 200
151 1 201 202
152 2 202 203
153 3 203 204
154 1 205 206
155 2 206 207
156 3 207 208
157 1 209 210
158 2 210 211
159 3 211 212
160 1 213 214
161 2 214 215
162 3 215 216
163 1 217 218
164 2 218 219
165 3 219 220
166 1 221 222
167 2 222 223
168 3 223 224
169 1 225 226
170 2 226 227
171 3 227 228
172 1 229 230
173 2 230 231
174 3 231 232
175 1 233 234
176 2 234 235
177 3 235 236
178 1 237 238
179 2 238 239
180 3 239 240
181 1 241 242
182 2 242 243
183 3 243 244
184 1 245 246
185 2 246 247
186 3 247 248
187 1 249 250
188 2 250 251
189 3 251 252
190 1 253 254
191 2 254 255
192 3 255 256
193 1 257 258
194 2 258 259
195 3 259 260
196 1 261 262
197 2 262 263
198 3 263 264
199 1 265 266
200 2 266 267
201 3 267 268
202 1 269 270
203 2 270 271
204 3 271 272
205 1 273 274
206 2 274 275
207 3 275 276
208 1 277 278
209 2 278 279
210 3 279 280
211 1 281 282
212 2 282 283
213 3 283 284
214 1 285 286
215 2 286 287
216 3 287 288
Angles
1 1 1 2 3
2 1 2 3 4
3 1 5 6 7
4 1 6 7 8
5 1 9 10 11
6 1 10 11 12
7 1 13 14 15
8 1 14 15 16
9 1 17 18 19
10 1 18 19 20
11 1 21 22 23
12 1 22 23 24
13 1 25 26 27
14 1 26 27 28
15 1 29 30 31
16 1 30 31 32
17 1 33 34 35
18 1 34 35 36
19 1 37 38 39
20 1 38 39 40
21 1 41 42 43
22 1 42 43 44
23 1 45 46 47
24 1 46 47 48
25 1 49 50 51
26 1 50 51 52
27 1 53 54 55
28 1 54 55 56
29 1 57 58 59
30 1 58 59 60
31 1 61 62 63
32 1 62 63 64
33 1 65 66 67
34 1 66 67 68
35 1 69 70 71
36 1 70 71 72
37 1 73 74 75
38 1 74 75 76
39 1 77 78 79
40 1 78 79 80
41 1 81 82 83
42 1 82 83 84
43 1 85 86 87
44 1 86 87 88
45 1 89 90 91
46 1 90 91 92
47 1 93 94 95
48 1 94 95 96
49 1 97 98 99
50 1 98 99 100
51 1 101 102 103
52 1 102 103 104
53 1 105 106 107
54 1 106 107 108
55 1 109 110 111
56 1 110 111 112
57 1 113 114 115
58 1 114 115 116
59 1 117 118 119
60 1 118 119 120
61 1 121 122 123
62 1 122 123 124
63 1 125 126 127
64 1 126 127 128
65 1 129 130 131
66 1 130 131 132
67 1 133 134 135
68 1 134 135 136
69 1 137 138 139
70 1 138 139 140
71 1 141 142 143
72 1 142 143 144
73 1 145 146 147
74 1 146 147 148
75 1 149 150 151
76 1 150 151 152
77 1 153 154 155
78 1 154 155 156
79 1 157 158 159
80 1 158 159 160
81 1 161 162 163
82 1 162 163 164
83 1 165 166 167
84 1 166 167 168
85 1 169 170 171
86 1 170 171 172
87 1 173 174 175
88 1 174 175 176
89 1 177 178 179
90 1 178 179 180
91 1 181 182 183
92 1 182 183 184
93 1 185 186 187
94 1 186 187 188
95 1 189 190 191
96 1 190 191 192
97 1 193 194 195
98 1 194 195 196
99 1 197 198 199
100 1 198 199 200
101 1 201 202 203
102 1 202 203 204
103 1 205 206 207
104 1 206 207 208
105 1 209 210 211
106 1 210 211 212
107 1 213 214 215
108 1 214 215 216
109 1 217 218 219
110 1 218 219 220
111 1 221 222 223
112 1 222 223 224
113 1 225 226 227
114 1 226 227 228
115 1 229 230 231
116 1 230 231 232
117 1 233 234 235
118 1 234 235 236
119 1 237 238 239
120 1 238 239 240
121 1 241 242 243
122 1 242 243 244
123 1 245 246 247
124 1 246 247 248
125 1 249 250 251
126 1 250 251 252
127 1 253 254 255
128 1 254 255 256
129 1 257 258 259
130 1 258 259 260
131 1 261 262 263
132 1 262 263 264
133 1 265 266 267
134 1 266 267 268
135 1 269 270 271
136 1 270 271 272
137 1 273 274 275
138 1 274 275 276
139 1 277 278 279
140 1 278 279 280
141 1 281 282 283
142 1 282 283 284
143 1 285 286 287
144 1 286 287 288

View File

@ -0,0 +1,75 @@
# Simple virus model.
#
units lj
dimension 3
boundary f f f
atom_style angle
pair_style lj/cut 4.0
pair_modify shift yes
bond_style harmonic
angle_style harmonic
neigh_modify binsize 5.0
read_data init.data
group NC type 1
group rest type 2 3 4
group CA type 3
timestep 0.0001
variable R equal "ramp(v_R0,v_R1)"
variable RF equal ${R1}
variable Rw equal "v_RF - 0.5"
print "**** R scales from ${R0} to ${R1} and wall is at ${Rw} ****"
variable rr atom "sqrt( x*x + y*y + z*z )"
compute ravg CA reduce ave v_rr
region my_sphere sphere 0 0 0 v_Rw side out
fix wall all wall/region my_sphere lj126 1.0 1.0 1.1225
variable U equal pe
variable epair equal epair
fix U all ave/time 2000 1 2000 v_U v_R v_epair
variable t equal time
variable ravg equal c_ravg
variable T equal 0.25
variable pU equal f_U[1]
variable pE equal f_U[3]
fix out all print 2000 "$t ${pU} $R ${ravg} ${pE}" file thermo.${DATA}.dat screen no
fix step1 rest nvt temp $T $T 1.0
fix step2 NC nvt/manifold/rattle 1e-4 10 sphere v_R temp $T $T 1.0
special_bonds lj 0 0 0
angle_coeff 1 250.0 180.0
thermo_style custom time step pe v_pU ke etotal epair temp v_R
thermo 10000
compute pe all pe/atom
fix pe_avg all ave/atom 100 50 5000 c_pe
dump traj_all all custom 5000 virus.dump id type x y z f_pe_avg
run 2500000
unfix out
fix step1 rest nvt temp $T $T 1.0
fix step2 NC nvt/manifold/rattle 1e-4 10 sphere ${R1} temp $T $T 1.0
run 1000000

View File

@ -0,0 +1,35 @@
# Install/unInstall package files in LAMMPS
# mode = 0/1/2 for uninstall/install/update
#
# Copied from USER-MISC.
#
mode=$1
# arg1 = file, arg2 = file it depends on
action () {
if (test $mode = 0) then
rm -f ../$1
elif (! cmp -s $1 ../$1) then
if (test -z "$2" || test -e ../$2) then
cp $1 ..
if (test $mode = 2) then
echo " updating src/$1"
fi
fi
elif (test ! -n "$2") then
if (test ! -e ../$2) then
rm -f ../$1
fi
fi
}
# all package files
# Factory needs to be remade after each modification/new manifold,
# other than that it is all good.
for file in *.cpp *.h; do
action $file
# manifold_factory.cpp
done

89
src/USER-MANIFOLD/README Normal file
View File

@ -0,0 +1,89 @@
+==============================================================================+
This file is a part of the USER-MANIFOLD package.
This package allows LAMMPS to perform MD simulations of particles
constrained on a manifold (i.e., a 2D subspace of the 3D simulation
box). It achieves this using the RATTLE constraint algorithm applied
to single-particle constraint functions g(xi,yi,zi) = 0 and their
derivative (i.e. the normal of the manifold) n = grad(g).
Stefan Paquay, s.paquay@tue.nl
Applied Physics/Theory of Polymers and Soft Matter,
Eindhoven University of Technology (TU/e), The Netherlands
Thanks to Remy Kusters at TU/e for testing.
This software is distributed under the GNU General Public License.
+==============================================================================+
At the moment we have a few manifolds available, extending them is very easy:
To add a new manifold, do the following in the "USER-MANIFOLD" directory:
0. Create a new pair of source/header files, and name them "manifold_*.cpp/h",
where you should replace '*' with some (descriptive) name.
1. In the header file, add the following:
a. Include guards.
b. Make sure you #include "manifold.h"
c. In namespace LAMMPS_NS, add a new class that inherits publicly from manifold
and protectedly from Pointers.
d. The constructor has to take ( LAMMPS*, int, char ** ) as arguments,
and should initialise Pointers.
e. The header file has to contain somewhere the macro ManifoldStyle with as
first argument the name of the manifold and as second argument the name
of the class implementing this manifold. The macro expands into some code
that registers the manifold during static initialisation, before main is
entered.
2. In the source file, make sure you implement the following (of course,
you can also implement these in the header):
+====================================+=========================================+
| Function signature | Purpose |
+====================================+=========================================+
| destructor | Free space (can be empty) |
| constructor | Has to take (LAMMPS *lmp, int, char **) |
| | as arguments and has to initialise |
| | Pointers with lmp. |
| double g( double *x ) | A function that is 0 if the 3-vector |
| | x is on the manifold. |
| void n( double *x, double *n ) | Stores the normal of position x in the |
| | 3-vector n. |
| static const char *type() | Returns text that identifies the |
| | manifold in LAMMPS input scripts. |
| const char *id() | Should return whatever type() returns. |
| static int expected_argc() | Returns the number of arguments needed |
| | for the construction/initialisation of |
| | your manifold. Example: Sphere only |
| | needs a radius, so it returns 1. The |
| | spine needs 5 parameters, so it |
| | returns 5. |
| int nparams() | Should return same as expected_argc() |
+====================================+=========================================+
If the above instructions seem a bit intimidating, you can get by just fine
by copying an existing manifold and modifying it. See e.g. manifold_sphere for
a relatively simple manifold.
With those things in place, the install script should be able to add your
manifold to LAMMPS without any extra work, so just running
make yes-user-manifold
make <your_architecture>
should (re)compile LAMMPS with the manifolds added.
+==============================================================================+
Obviously, you need some parameters to represent the surface, such as the radius
in the case of a sphere. These should be passed to the nve/manifold/rattle fix,
with the following syntax:
fix ID group-ID nve/manifold/rattle tol maxit manifold_style args
tol = tolerance to which RATTLE tries to satisfy the constraints
maxit = maximum number of iterations RATTLE uses each time step
manifold_style = manifold style, should equal what type() of the desired
manifold returns
args = parameters for the manifold, order is manifold-dependent.
can be equal style variables
+==============================================================================+

View File

@ -0,0 +1,184 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "math.h"
#include "string.h"
#include "stdlib.h"
#include "atom.h"
#include "update.h"
#include "respa.h"
#include "error.h"
#include "force.h"
#include "manifold.h"
#include "fix_manifoldforce.h" // For stuff
#include "manifold_factory.h" // For constructing manifold
using namespace LAMMPS_NS;
using namespace FixConst;
// Helper functions for parameters/equal style variables in input script
inline bool was_var( const char *arg )
{
return strstr( arg, "v_" ) == arg;
}
inline bool str_eq( const char *str1, const char *str2 )
{
return strcmp(str1,str2) == 0;
}
/* ---------------------------------------------------------------------- */
FixManifoldForce::FixManifoldForce(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
int me = -1;
MPI_Comm_rank(world,&me);
// Check the min-style:
int good_minner = str_eq(update->minimize_style,"hftn") |
str_eq(update->minimize_style,"quickmin");
if( !good_minner){
error->warning(FLERR,"Minimizing with fix manifoldforce without hftn or quickmin is fishy");
}
// Command is given as
// fix <name> <group> manifoldforce manifold_name manifold_args
if( narg < 5 ){
error->all(FLERR,"Illegal fix manifoldforce! No manifold given");
}
const char *m_name = arg[3];
ptr_m = create_manifold(m_name,lmp,narg,arg);
// Construct manifold from factory:
if( !ptr_m ){
char msg[2048];
snprintf(msg,2048,"Manifold pointer for manifold '%s' was NULL for some reason", arg[3]);
error->all(FLERR,msg);
}
// After constructing the manifold, you can safely make
// room for the parameters
nvars = ptr_m->nparams();
if( narg < nvars+4 ){
char msg[2048];
sprintf(msg,"Manifold %s needs at least %d argument(s)!",
m_name, nvars);
error->all(FLERR,msg);
}
*(ptr_m->get_params()) = new double[nvars];
if( ptr_m->get_params() == NULL ){
error->all(FLERR,"Parameter pointer was NULL!");
}
// This part here stores the names/text of each argument,
// determines which params are equal-style variables,
// and sets the values of those arguments that were _not_
// equal style vars (so that they are not overwritten each time step).
double *params = *(ptr_m->get_params());
for( int i = 0; i < nvars; ++i ){
if( was_var( arg[i+4] ) )
error->all(FLERR,"Equal-style variables not allowed with fix manifoldforce");
// Use force->numeric to trigger an error if arg is not a number.
params[i] = force->numeric(FLERR,arg[i+4]);
}
// Perform any further initialisation for the manifold that depends on params:
ptr_m->post_param_init();
}
/* ---------------------------------------------------------------------- */
int FixManifoldForce::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= MIN_POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixManifoldForce::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
int nlevels_respa = ((Respa *) update->integrate)->nlevels;
for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {
((Respa *) update->integrate)->copy_flevel_f(ilevel);
post_force_respa(vflag,ilevel,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel);
}
}
}
/* ---------------------------------------------------------------------- */
void FixManifoldForce::min_setup(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixManifoldForce::post_force(int vflag)
{
double **x = atom->x;
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double n[3];
double invn2;
double dot;
for (int i = 0; i < nlocal; i++){
if (mask[i] & groupbit) {
// Determine normal of particle:
ptr_m->n(x[i],n);
invn2 = 1.0 / ( n[0]*n[0] + n[1]*n[1] + n[2]*n[2] );
dot = f[i][0]*n[0] + f[i][1]*n[1] + f[i][2]*n[2];
f[i][0] -= dot*n[0] * invn2;
f[i][1] -= dot*n[1] * invn2;
f[i][2] -= dot*n[2] * invn2;
}
}
}
/* ---------------------------------------------------------------------- */
void FixManifoldForce::post_force_respa(int vflag, int ilevel, int iloop)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixManifoldForce::min_post_force(int vflag)
{
post_force(vflag);
}

View File

@ -0,0 +1,81 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------
This file is a part of the USER-MANIFOLD package.
Copyright (2013-2015) Stefan Paquay, Eindhoven University of Technology.
License: GNU General Public License.
See the README file in the top-level LAMMPS directory.
This file is part of the user-manifold package written by
Stefan Paquay at the Eindhoven University of Technology.
This module makes it possible to do MD with particles constrained
to pretty arbitrary manifolds characterised by some constraint function
g(x,y,z) = 0 and its normal grad(g). The number of manifolds available
right now is limited but can be extended straightforwardly by making
a new class that inherits from manifold and implements all pure virtual
methods.
This fix subtracts force components that point out of the manifold,
useful for minimisation on surfaces.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(manifoldforce,FixManifoldForce)
#else
#ifndef LMP_FIX_MANIFOLDFORCE_H
#define LMP_FIX_MANIFOLDFORCE_H
#include "fix.h"
#include "manifold.h"
namespace LAMMPS_NS {
class FixManifoldForce : public Fix {
public:
FixManifoldForce(class LAMMPS *, int, char **);
int setmask();
void setup(int);
void min_setup(int);
void post_force(int);
void post_force_respa(int, int, int);
void min_post_force(int);
private:
manifold *ptr_m;
// Stuff to store the parameters in.
int nvars; // # of args after manifold name.
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
*/

View File

@ -0,0 +1,640 @@
/* ----------------------------------------------------------------------
Lammps - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
-----------------------------------------------------------------------
This file is a part of the USER-MANIFOLD package.
Copyright (2013-2014) Stefan Paquay, Eindhoven University of Technology.
License: GNU General Public License.
See the README file in the top-level LAMMPS directory.
This file is part of the user-manifold package written by
Stefan Paquay at the Eindhoven University of Technology.
This module makes it possible to do MD with particles constrained
to pretty arbitrary manifolds characterised by some constraint function
g(x,y,z) = 0 and its normal grad(g). The number of manifolds available
right now is limited but can be extended straightforwardly by making
a new class that inherits from manifold and implements all pure virtual
methods.
Thanks to Remy Kusters for beta-testing!
------------------------------------------------------------------------- */
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "atom.h"
#include "force.h"
#include "update.h"
#include "respa.h"
#include "error.h"
#include "group.h"
#include "math.h"
#include "input.h"
#include "variable.h"
#include "citeme.h"
#include "memory.h"
#include "comm.h"
#include "fix_nve_manifold_rattle.h"
#include "manifold_factory.h"
#include "manifold.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum { CONST, EQUAL }; // For treating the variables.
static const char* cite_fix_nve_manifold_rattle =
"fix nve/manifold/rattle command:\n\n"
"@article{paquay-2016,\n"
" author = {Paquay, Stefan and Kusters, Remy}, \n"
" eprint = {arXiv:1411.3019}, \n"
" title = {A method for molecular dynamics on curved surfaces}, \n"
" url = {http://arxiv.org/abs/1411.3019}, \n"
" year = 2016, \n"
" journal = {Arxiv preprint, \\url{http://arxiv.org/abs/1411.3019}}, \n"
" note = {To be published in Biophys. J.}, \n"
"}\n\n";
/* -----------------------------------------------------------------------------
---------------------------------------------------------------------------*/
FixNVEManifoldRattle::FixNVEManifoldRattle( LAMMPS *lmp, int &narg, char **arg,
int error_on_unknown_keyword )
: Fix(lmp,narg,arg)
{
if( lmp->citeme) lmp->citeme->add(cite_fix_nve_manifold_rattle);
if( narg < 6 ) error->all(FLERR, "Illegal fix nve/manifold/rattle command");
// Set all bits/settings:
time_integrate = 1;
dynamic_group_allow = 1;
size_vector = 0;
dof_flag = 1;
nevery = 0;
dtv = dtf = 0;
tolerance = force->numeric( FLERR, arg[3] );
max_iter = force->numeric( FLERR, arg[4] );
ptr_m = create_manifold(arg[5], lmp, narg, arg);
if( !ptr_m ){
error->all(FLERR,"Error creating manifold pointer");
}
nvars = ptr_m->nparams();
tstrs = new char*[nvars];
tvars = new int[nvars];
tstyle = new int[nvars];
is_var = new int[nvars];
if( !tstrs || !tvars || !tstyle || !is_var ){
error->all(FLERR, "Error creating manifold arg arrays");
}
// Loop over manifold args:
for( int i = 0; i < nvars; ++i ){
int len = 0, offset = 0;
if( was_var( arg[i+6] ) ){
len = strlen(arg[i+6]) - 1; // -1 because -2 for v_, +1 for \0.
is_var[i] = 1;
offset = 2;
}else{
force->numeric(FLERR,arg[i+6]); // Check if legal number.
len = strlen( arg[i+6] ) + 1; // +1 for \0.
is_var[i] = 0;
}
tstrs[i] = new char[len];
if( tstrs[i] == NULL ) error->all(FLERR,"Error allocating space for args.");
strcpy( tstrs[i], arg[i+6] + offset );
}
*ptr_m->get_params() = new double[nvars];
if( !(*ptr_m->get_params()) ) error->all(FLERR,"Failed to allocate params!");
for( int i = 0; i < nvars; ++i ){
// If param i was variable type, it will be set later...
(*ptr_m->get_params())[i] = is_var[i] ? 0.0 : force->numeric( FLERR, arg[i+6] );
}
ptr_m->post_param_init();
// Loop over rest of args:
int argi = 6 + nvars;
while( argi < narg ){
if( strcmp(arg[argi], "every") == 0 ){
nevery = force->inumeric(FLERR,arg[argi+1]);
argi += 2;
}else if( error_on_unknown_keyword ){
char msg[2048];
sprintf(msg,"Error parsing arg \"%s\".\n", arg[argi]);
error->all(FLERR, msg);
}else{
argi += 1;
}
}
}
/* -----------------------------------------------------------------------------
---------------------------------------------------------------------------*/
FixNVEManifoldRattle::~FixNVEManifoldRattle()
{
if( tstrs ){
for( int i = 0; i < nvars; ++i ){
delete [] tstrs[i];
}
delete [] tstrs;
}
if( tvars ) delete [] tvars;
if( tstyle ) delete [] tstyle;
if( is_var ) delete [] is_var;
}
/* -----------------------------------------------------------------------------
---------------------------------------------------------------------------*/
void FixNVEManifoldRattle::reset_dt()
{
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
}
void FixNVEManifoldRattle::print_stats( const char *header )
{
double n = stats.natoms;
if( n > 0 ){
stats.x_iters_per_atom += stats.x_iters / n;
stats.v_iters_per_atom += stats.v_iters / n;
}
double x_iters = 0, v_iters = 0;
bigint ntimestep = update->ntimestep;
int me = -1;
MPI_Comm_rank(world,&me);
MPI_Allreduce(&stats.x_iters_per_atom,&x_iters,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&stats.v_iters_per_atom,&v_iters,1,MPI_DOUBLE,MPI_SUM,world);
// Set iters back to zero:
stats.x_iters_per_atom = stats.x_iters = 0;
stats.v_iters_per_atom = stats.v_iters = 0;
if( me == 0 ){
double inv_tdiff = 1.0/( static_cast<double>(ntimestep) - stats.last_out );
stats.last_out = ntimestep;
fprintf(screen, "%s stats for time step " BIGINT_FORMAT " on %d atoms:\n",
header, ntimestep, stats.natoms);
fprintf(screen, " iters/atom: x = %f, v = %f, dofs removed %d",
x_iters * inv_tdiff, v_iters * inv_tdiff, stats.dofs_removed);
fprintf(screen,"\n");
}
}
/* -----------------------------------------------------------------------------
---------------------------------------------------------------------------*/
int FixNVEManifoldRattle::was_var( const char *str )
{
if( strlen(str) > 2 ){
return (str[0] == 'v') && (str[1] == '_');
}else{
return 0;
}
}
/* -----------------------------------------------------------------------------
---------------------------------------------------------------------------*/
int FixNVEManifoldRattle::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
if( nevery > 0 ) mask |= END_OF_STEP;
return mask;
}
/* -----------------------------------------------------------------------------
---------------------------------------------------------------------------*/
void FixNVEManifoldRattle::init()
{
// Makes sure the manifold params are set initially.
update_var_params();
reset_dt();
}
void FixNVEManifoldRattle::update_var_params()
{
if( nevery > 0 ){
stats.x_iters = 0;
stats.v_iters = 0;
stats.natoms = 0;
stats.x_iters_per_atom = 0.0;
stats.v_iters_per_atom = 0.0;
}
double **ptr_params = ptr_m->get_params();
for( int i = 0; i < nvars; ++i ){
if( is_var[i] ){
tvars[i] = input->variable->find(tstrs[i]);
if( tvars[i] < 0 ){
error->all(FLERR,
"Variable name for fix nve/manifold/rattle does not exist");
}
if( input->variable->equalstyle(tvars[i]) ){
tstyle[i] = EQUAL;
double new_val = input->variable->compute_equal(tvars[i]);
// fprintf( stdout, "New value of var %d is now %f\n", i+1, new_val );
*(ptr_params[i]) = new_val;
}else{
error->all(FLERR,
"Variable for fix nve/manifold/rattle is invalid style");
}
}
}
}
/* -----------------------------------------------------------------------------
---------------------------------------------------------------------------*/
int FixNVEManifoldRattle::dof(int igroup)
{
int *mask = atom->mask;
int nlocal = atom->nlocal;
int natoms = 0;
for( int i = 0; i < nlocal; ++i ){
if(mask[i] & groupbit) ++natoms;
}
int dofs;
MPI_Allreduce( &natoms, &dofs, 1, MPI_INT, MPI_SUM, world );
// Make sure that, if there is just no or one atom, no dofs are subtracted,
// since for the first atom already 3 dofs are subtracted because of the
// centre of mass corrections:
if( dofs <= 1 ) dofs = 0;
stats.dofs_removed = dofs;
return dofs;
}
/* -----------------------------------------------------------------------------
---------------------------------------------------------------------------*/
double FixNVEManifoldRattle::memory_usage()
{
double bytes = 0.0;
bytes += sizeof(statistics);
bytes += sizeof(*ptr_m) + sizeof(ptr_m);
bytes += nvars*sizeof(double) + sizeof(double*);
bytes += nvars*( sizeof(char*) + 3*sizeof(int) );
return bytes;
}
/* -----------------------------------------------------------------------------
---------------------------------------------------------------------------*/
void FixNVEManifoldRattle::initial_integrate(int vflag)
{
update_var_params();
nve_x_rattle(igroup, groupbit);
}
/* -----------------------------------------------------------------------------
---------------------------------------------------------------------------*/
void FixNVEManifoldRattle::final_integrate()
{
nve_v_rattle(igroup, groupbit);
}
/* -----------------------------------------------------------------------------
---------------------------------------------------------------------------*/
void FixNVEManifoldRattle::end_of_step()
{
print_stats( "nve/manifold/rattle" );
}
/* -----------------------------------------------------------------------------
---------------------------------------------------------------------------*/
void FixNVEManifoldRattle::nve_x_rattle(int igroup, int groupbit)
{
double dtfm;
// update v and x of atoms in group
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int natoms = 0;
if (igroup == atom->firstgroup){
nlocal = atom->nfirst;
}
if (rmass) {
for (int i = 0; i < nlocal; i++){
if (mask[i] & groupbit){
natoms++;
dtfm = dtf / rmass[i];
rattle_manifold_x( x[i], v[i], f[i], dtv, dtfm, atom->tag[i] );
}
}
} else {
for (int i = 0; i < nlocal; i++){
if (mask[i] & groupbit) {
natoms++;
dtfm = dtf / mass[type[i]];
rattle_manifold_x( x[i], v[i], f[i], dtv, dtfm, atom->tag[i] );
}
}
}
if( nevery > 0 ){
// Count ALL atoms this fix works on:
MPI_Allreduce(&natoms,&stats.natoms,1,MPI_INT,MPI_SUM,world);
}
}
/* -----------------------------------------------------------------------------
---------------------------------------------------------------------------*/
void FixNVEManifoldRattle::nve_v_rattle(int igroup, int groupbit)
{
double dtfm;
// update v of atoms in group
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
rattle_manifold_v( v[i], f[i], x[i], dtfm, atom->tag[i] );
}
}
} else {
for (int i = 0; i < nlocal; i++){
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
rattle_manifold_v( v[i], f[i], x[i], dtfm, atom->tag[i] );
}
}
}
}
/* -----------------------------------------------------------------------------
---------------------------------------------------------------------------*/
void FixNVEManifoldRattle::rattle_manifold_x(double *x, double *v,
double *f, double dtv,
double dtfm, tagint tagi )
{
/*
A RATTLE update for the position constraint.
Original update is x += dtv * v_1/2
Now you do
v_1/2(lambda) = v_0 + dtfm * ( f + lambda*n_old )
and solve
xold - xnew + dtv * v_1/2(lambda) = 0
g(xnew) = 0
for x and lambda. The lambda you find then gives v_1/2 as well.
*/
double xo[3]; // Previous position to update from.
double vo[3]; // Previous velocity to update from.
double l = 0; // Lagrangian multiplier for constraint forces.
double R[4]; // System that is 0.
double dx[4]; // Update that follows from Newton iteration.
double no[3]; // Normal at xo.
double nn[3]; // Normal at x, the new position.
double res; // Residual.
int iters = 0; // Iterations used
double c = dtfm*dtv; // Used for iterating in the Newton loop:
double no_nn, nn_R;
vo[0] = v[0];
vo[1] = v[1];
vo[2] = v[2];
xo[0] = x[0];
xo[1] = x[1];
xo[2] = x[2];
double gg = ptr_m->g_and_n(x,no);
nn[0] = no[0];
nn[1] = no[1];
nn[2] = no[2];
double vt[3];
vt[0] = vo[0] + dtfm*f[0];
vt[1] = vo[1] + dtfm*f[1];
vt[2] = vo[2] + dtfm*f[2];
double no_dt[3];
no_dt[0] = dtfm*no[0];
no_dt[1] = dtfm*no[1];
no_dt[2] = dtfm*no[2];
// Assume that no_nn is roughly constant during iteration:
const double c_inv = 1.0 / c;
while ( 1 ) {
v[0] = vt[0] - l*no_dt[0];
v[1] = vt[1] - l*no_dt[1];
v[2] = vt[2] - l*no_dt[2];
R[0] = xo[0] - x[0] + dtv * v[0];
R[1] = xo[1] - x[1] + dtv * v[1];
R[2] = xo[2] - x[2] + dtv * v[2];
R[3] = gg;
// Analytic solution to system J*(dx,dy,dz,dl)^T = R
// no_nn = no[0]*nn[0] + no[1]*nn[1] + no[2]*nn[2];
nn_R = nn[0]*R[0] + nn[1]*R[1] + nn[2]*R[2];
no_nn = no[0]*nn[0] + no[1]*nn[1] + no[2]*nn[2];
double n_inv = 1.0 / no_nn;
// fprintf( screen, "nn_R = %f, no_nn = %f\n", nn_R, no_nn );
dx[3] = -nn_R - R[3];
dx[3] *= n_inv;
dx[0] = -R[0] - no[0]*dx[3];
dx[1] = -R[1] - no[1]*dx[3];
dx[2] = -R[2] - no[2]*dx[3];
dx[3] *= c_inv;
x[0] -= dx[0];
x[1] -= dx[1];
x[2] -= dx[2];
l -= dx[3];
res = infnorm<4>(R);
++iters;
if( (res < tolerance) || (iters >= max_iter) ) break;
// Update nn and g.
gg = ptr_m->g(x);
ptr_m->n(x,nn);
// gg = ptr_m->g(x);
}
if( iters >= max_iter && res > tolerance ){
char msg[2048];
sprintf(msg,"Failed to constrain atom %d (x = (%f, %f, %f)! res = %e, iters = %d\n",
tagi, x[0], x[1], x[2], res, iters);
error->one(FLERR,msg);
}
// "sync" x and v:
v[0] = vt[0] - l*no_dt[0];
v[1] = vt[1] - l*no_dt[1];
v[2] = vt[2] - l*no_dt[2];
stats.x_iters += iters;
}
/* -----------------------------------------------------------------------------
---------------------------------------------------------------------------*/
void FixNVEManifoldRattle::rattle_manifold_v(double *v, double *f,
double *x, double dtfm,
tagint tagi )
{
/*
The original update was
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
Now you add the rattle-like update:
vold - vnew + dtfm * F + mu * n_new = 0
dot( vnew, n_new ) = 0
*/
double vo[3]; // V at t + 1/2 dt
double l = 0; // Lagrangian multiplier for constraint forces.
double R[4]; // System that is 0.
double dv[4]; // Update that follows from Newton iteration.
double n[3]; // Normal.
double res; // Residual.
int iters = 0; // Iterations used
double c = dtfm; // Used for iterating in the Newton loop:
double nn2, nn_R;
vo[0] = v[0];
vo[1] = v[1];
vo[2] = v[2];
// Initial guess is unconstrained update:
v[0] += dtfm*f[0];
v[1] += dtfm*f[1];
v[2] += dtfm*f[2];
ptr_m->n(x,n);
double vt[3];
vt[0] = vo[0] + dtfm*f[0];
vt[1] = vo[1] + dtfm*f[1];
vt[2] = vo[2] + dtfm*f[2];
double no_dt[3];
no_dt[0] = dtfm*n[0];
no_dt[1] = dtfm*n[1];
no_dt[2] = dtfm*n[2];
nn2 = n[0]*n[0] + n[1]*n[1] + n[2]*n[2];
const double n_inv = 1.0 / nn2;
const double c_inv = 1.0 / c;
do{
R[0] = vt[0] - v[0] - l * no_dt[0];
R[1] = vt[1] - v[1] - l * no_dt[1];
R[2] = vt[2] - v[2] - l * no_dt[2];
R[3] = v[0]*n[0] + v[1]*n[1] + v[2]*n[2];
// Analytic solution to system J*(dx,dy,dz,dl)^T = R
nn_R = n[0]*R[0] + n[1]*R[1] + n[2]*R[2];
dv[3] = -nn_R - R[3];
dv[3] *= n_inv;
dv[0] = -n[0]*dv[3] - R[0];
dv[1] = -n[1]*dv[3] - R[1];
dv[2] = -n[2]*dv[3] - R[2];
dv[3] *= c_inv;
v[0] -= dv[0];
v[1] -= dv[1];
v[2] -= dv[2];
l -= dv[3];
res = infnorm<4>(R);
++iters;
}while( (res > tolerance) && (iters < max_iter) );
if( iters >= max_iter && res >= tolerance ){
char msg[2048];
sprintf(msg,"Failed to constrain atom %d (x = (%f, %f, %f)! res = %e, iters = %d\n",
tagi, x[0], x[1], x[2], res, iters);
error->all(FLERR,msg);
}
stats.v_iters += iters;
}

View File

@ -0,0 +1,170 @@
/* ----------------------------------------------------------------------
Lammps - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
-----------------------------------------------------------------------
This file is a part of the USER-MANIFOLD package.
Copyright (2013-2014) Stefan Paquay, Eindhoven University of Technology.
License: GNU General Public License.
See the README file in the top-level LAMMPS directory.
This file is part of the user-manifold package written by
Stefan Paquay at the Eindhoven University of Technology.
This module makes it possible to do MD with particles constrained
to pretty arbitrary manifolds characterised by some constraint function
g(x,y,z) = 0 and its normal grad(g). The number of manifolds available
right now is limited but can be extended straightforwardly by making
a new class that inherits from manifold and implements all pure virtual
methods.
Thanks to Remy Kusters for beta-testing!
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(nve/manifold/rattle,FixNVEManifoldRattle)
#else
#ifndef LMP_FIX_NVE_MANIFOLD_RATTLE_H
#define LMP_FIX_NVE_MANIFOLD_RATTLE_H
#include "fix.h"
#include "manifold.h"
namespace LAMMPS_NS {
class FixNVEManifoldRattle : public Fix {
public:
struct statistics {
statistics() : x_iters(0), v_iters(0), x_iters_per_atom(0),
v_iters_per_atom(0), natoms(0), dofs_removed(0),
last_out(0) {}
double x_iters, v_iters;
double x_iters_per_atom;
double v_iters_per_atom;
int natoms;
int dofs_removed;
bigint last_out;
};
FixNVEManifoldRattle(LAMMPS *, int &, char **, int = 1);
virtual ~FixNVEManifoldRattle();
// All this stuff is interface, so you DO need to implement them.
// Just delegate them to the workhorse classes.
virtual int setmask();
virtual void initial_integrate(int);
virtual void final_integrate();
virtual void init();
virtual void reset_dt();
virtual void end_of_step();
virtual int dof(int);
virtual void setup(int){} // Not needed for fixNVE but is for fixNVT
virtual double memory_usage();
protected:
int nevery;
double dtv, dtf;
double tolerance;
int max_iter;
char **tstrs;
int *tvars;
int *tstyle;
int *is_var;
statistics stats;
int update_style;
int nvars;
manifold *ptr_m;
void print_stats( const char * );
int was_var( const char * );
virtual void update_var_params();
virtual void rattle_manifold_x( double *, double *, double *, double, double, tagint );
virtual void rattle_manifold_v( double *, double *, double *, double, tagint );
virtual void nve_x_rattle(int, int);
virtual void nve_v_rattle(int, int);
};
}
#endif // LMP_FIX_NVE_MANIFOLD_RATTLE_H
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: There is no manifold named ...
Self-explanatory. You requested a manifold whose name was not
registered at the factory.
E: Manifold pointer was NULL for some reason!
This indicates a bug. The factory was unable to properly create
the requested manifold even though it was registered. Send the
maintainer an e-mail.
E: Manifold ... needs at least ... argument(s)!
Self-explanatory. Provide enough arguments for the proper
creating of the requested manifold.
E: Parameter pointer was NULL!
This indicates a bug. The array that contains the parameters
could not be allocated. Send the maintainer an e-mail.
E: Could not allocate space for arg!
One of the arguments provided was too long (it contained
too many characters)
E: Option ... needs ... argument(s):
Self-explanatory. Read the documentation of this fix properly.
E: Illegal fix nve/manifold/rattle command! Option ... not recognized!
Self-explanatory. The given option(s) do not exist.
E: Variable name for fix nve/manifold/rattle does not exist
Self-explanatory.
E: Variable for fix nve/manifold/rattle is invalid style
fix nve/manifold/rattle only supports equal style variables.
*/

View File

@ -0,0 +1,415 @@
/* ----------------------------------------------------------------------
Lammps - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
-----------------------------------------------------------------------
This file is a part of the USER-MANIFOLD package.
Copyright (2013-2014) Stefan Paquay, Eindhoven University of Technology.
License: GNU General Public License.
See the README file in the top-level LAMMPS directory.
This file is part of the user-manifold package written by
Stefan Paquay at the Eindhoven University of Technology.
This module makes it possible to do MD with particles constrained
to pretty arbitrary manifolds characterised by some constraint function
g(x,y,z) = 0 and its normal grad(g). The number of manifolds available
right now is limited but can be extended straightforwardly by making
a new class that inherits from manifold and implements all pure virtual
methods.
Thanks to Remy Kusters for beta-testing!
------------------------------------------------------------------------- */
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "atom.h"
#include "force.h"
#include "update.h"
#include "respa.h"
#include "error.h"
#include "group.h"
#include "math.h"
#include "input.h"
#include "variable.h"
#include "citeme.h"
#include "memory.h"
#include "comm.h"
#include "modify.h"
#include "compute.h"
#include "fix_nvt_manifold_rattle.h"
#include "manifold.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum {CONSTANT,EQUAL};
enum {NOBIAS,BIAS};
static const char* cite_fix_nvt_manifold_rattle =
"fix nve/manifold/rattle command:\n\n"
"@article{paquay-2016,\n"
" author = {Paquay, Stefan and Kusters, Remy}, \n"
" eprint = {arXiv:1411.3019}, \n"
" title = {A method for molecular dynamics on curved surfaces}, \n"
" url = {http://arxiv.org/abs/1411.3019}, \n"
" year = 2016, \n"
" journal = {Arxiv preprint, \\url{http://arxiv.org/abs/1411.3019}}, \n"
" note = {To be published in Biophys. J.}, \n"
"}\n\n";
/* ---------------------------------------------------------------------- */
FixNVTManifoldRattle::FixNVTManifoldRattle(LAMMPS *lmp, int narg, char **arg,
int error_on_unknown_keyword )
: FixNVEManifoldRattle(lmp,narg,arg, 0)
{
if (lmp->citeme) lmp->citeme->add(cite_fix_nvt_manifold_rattle);
if( narg < 6 ) error->all(FLERR,"Illegal fix nvt/manifold/rattle command");
// Set all bits/settings:
dof_flag = 1;
dthalf = dt4 = dt8 = 0;
t_start = t_stop = t_period = t_current = t_target = ke_target = 0.0;
t_freq = drag = tdrag_factor = 0;
boltz = force->boltz, nktv2p = force->nktv2p;
tdof = 0;
mtchain = 3;
factor_eta = 0.0;
which = got_temp = 0;
int argi = 6 + ptr_m->nparams();
while( argi < narg )
{
if( strcmp( arg[argi], "temp") == 0 ){
if( argi+3 >= narg )
error->all(FLERR,"Keyword 'temp' needs 3 arguments");
t_start = force->numeric(FLERR, arg[argi+1]);
t_stop = force->numeric(FLERR, arg[argi+2]);
t_period = force->numeric(FLERR, arg[argi+3]);
t_target = t_start;
got_temp = 1;
argi += 4;
}else if( strcmp( arg[argi], "tchain" ) == 0 ){
if( argi+1 >= narg )
error->all(FLERR,"Keyword 'tchain' needs 1 argument");
mtchain = force->inumeric(FLERR, arg[argi+1]);
argi += 2;
}else if( error_on_unknown_keyword ){
char msg[2048];
sprintf(msg,"Error parsing arg \"%s\".\n", arg[argi]);
error->all(FLERR, msg);
}else{
argi += 1;
}
}
reset_dt();
if( !got_temp ) error->all(FLERR,"Fix nvt/manifold/rattle needs 'temp'!");
if( t_period < 0.0 ){
error->all(FLERR,"Fix nvt/manifold/rattle damping parameter must be > 0.0");
}
// Create temperature compute:
const char *fix_id = arg[1];
int n = strlen(fix_id)+6;
id_temp = new char[n];
strcpy(id_temp,fix_id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = group->names[igroup];
newarg[2] = (char*) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
int icompute = modify->find_compute(id_temp);
if( icompute < 0 ){
error->all(FLERR,"Temperature ID for fix nvt/manifold/rattle "
"does not exist");
}
temperature = modify->compute[icompute];
if( temperature->tempbias ) which = BIAS;
else which = NOBIAS;
// Set t_freq from t_period
t_freq = 1.0 / t_period;
// Init Nosé-Hoover chain:
eta = new double[mtchain];
eta_dot = new double[mtchain+1];
eta_dotdot = new double[mtchain];
eta_mass = new double[mtchain];
eta_dot[mtchain] = 0.0;
eta_dot[mtchain] = 0.0;
for( int ich = 0; ich < mtchain; ++ich ){
eta[ich] = eta_dot[ich] = eta_dotdot[ich] = 0.0;
}
}
/* ---------------------------------------------------------------------- */
FixNVTManifoldRattle::~FixNVTManifoldRattle()
{
// Deallocate heap-allocated objects.
if( eta ) delete[] eta;
if( eta_dot ) delete[] eta_dot;
if( eta_dotdot ) delete[] eta_dotdot;
if( eta_mass ) delete[] eta_mass;
modify->delete_compute(id_temp);
if( id_temp ) delete[] id_temp;
}
int FixNVTManifoldRattle::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
if( nevery > 0 ) mask |= END_OF_STEP;
return mask;
}
/* --------------------------------------------------------------------------
Check that force modification happens before position and velocity update.
Make sure respa is not used.
------------------------------------------------------------------------- */
void FixNVTManifoldRattle::init()
{
// Makes sure the manifold params are set initially.
update_var_params();
int icompute = modify->find_compute(id_temp);
if( icompute < 0 ){
error->all(FLERR,"Temperature ID for fix nvt/manifold/rattle "
"does not exist");
}
temperature = modify->compute[icompute];
if( temperature->tempbias ) which = BIAS;
else which = NOBIAS;
}
void FixNVTManifoldRattle::setup(int vflag)
{
compute_temp_target();
t_current = temperature->compute_scalar();
tdof = temperature->dof;
// Compute/set eta-masses:
double inv_t_freq2 = 1.0 / (t_freq*t_freq);
eta_mass[0] = tdof * boltz * t_target * inv_t_freq2;
for( int ich = 1; ich < mtchain; ++ich ){
eta_mass[ich] = boltz * t_target * inv_t_freq2;
}
for( int ich = 1; ich < mtchain; ++ich ){
eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1] -
boltz * t_target ) / eta_mass[ich];
}
}
void FixNVTManifoldRattle::compute_temp_target()
{
t_current = temperature->compute_scalar();
tdof = temperature->dof;
double delta = update->ntimestep - update->beginstep;
if (delta != 0.0){
delta /= update->endstep - update->beginstep;
}
tdof = temperature->dof;
t_target = t_start + delta * (t_stop-t_start);
ke_target = tdof * boltz * t_target;
}
void FixNVTManifoldRattle::nhc_temp_integrate()
{
int ich;
// t_current = temperature->compute_scalar();
// tdof = temperature->dof;
compute_temp_target();
double expfac, kecurrent = tdof * boltz * t_current;
double inv_t_freq2 = 1.0 / (t_freq*t_freq);
eta_mass[0] = tdof * boltz * t_target * inv_t_freq2;
for( int ich = 1; ich < mtchain; ++ich ){
eta_mass[ich] = boltz * t_target * inv_t_freq2;
}
if( eta_mass[0] > 0.0 ){
eta_dotdot[0] = (kecurrent - ke_target)/eta_mass[0];
}else{
eta_dotdot[0] = 0;
}
for( ich = mtchain-1; ich > 0; --ich ){
expfac = exp(-dt8*eta_dot[ich+1]);
eta_dot[ich] *= expfac;
eta_dot[ich] += eta_dotdot[ich] * dt4;
eta_dot[ich] *= tdrag_factor * expfac;
}
expfac = exp(-dt8*eta_dot[1]);
eta_dot[0] *= expfac;
eta_dot[0] += eta_dotdot[0] * dt4;
eta_dot[0] *= tdrag_factor * expfac;
factor_eta = exp(-dthalf*eta_dot[0]);
if( factor_eta == 0 ){
char msg[2048];
sprintf(msg, "WTF, factor_eta is 0! dthalf = %f, eta_dot[0] = %f",
dthalf, eta_dot[0]);
error->all(FLERR,msg);
}
nh_v_temp();
t_current *= factor_eta*factor_eta;
kecurrent = tdof * boltz * t_current;
if( eta_mass[0] > 0.0 ){
eta_dotdot[0] = (kecurrent - ke_target) / eta_mass[0];
}else{
eta_dotdot[0] = 0.0;
}
for( int ich = 1; ich < mtchain; ++ich ){
eta[ich] += dthalf*eta_dot[ich];
}
eta_dot[0] *= expfac;
eta_dot[0] += eta_dotdot[0]*dt4;
eta_dot[0] *= expfac;
for( int ich = 1; ich < mtchain; ++ich ){
expfac = exp(-dt8*eta_dot[ich+1]);
eta_dot[ich] *= expfac;
eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1]
- boltz*t_target) / eta_mass[ich];
eta_dot[ich] *= eta_dotdot[ich] * dt4;
eta_dot[ich] *= expfac;
}
}
void FixNVTManifoldRattle::nh_v_temp()
{
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if( igroup == atom->firstgroup) nlocal = atom->nfirst;
if( which == NOBIAS ){
for( int i = 0; i < nlocal; ++i ){
if( mask[i] & groupbit ){
v[i][0] *= factor_eta;
v[i][1] *= factor_eta;
v[i][2] *= factor_eta;
}
}
}else if( which == BIAS ){
for( int i = 0; i < nlocal; ++i ){
if( mask[i] & groupbit ){
temperature->remove_bias(i,v[i]);
v[i][0] *= factor_eta;
v[i][1] *= factor_eta;
v[i][2] *= factor_eta;
temperature->restore_bias(i,v[i]);
}
}
}
}
// Most of this logic is based on fix_nh:
void FixNVTManifoldRattle::initial_integrate(int vflag)
{
update_var_params();
compute_temp_target();
nhc_temp_integrate();
nve_x_rattle(igroup, groupbit);
}
void FixNVTManifoldRattle::final_integrate()
{
nve_v_rattle(igroup, groupbit);
nhc_temp_integrate();
}
/* ---------------------------------------------------------------------- */
void FixNVTManifoldRattle::reset_dt()
{
FixNVEManifoldRattle::reset_dt();
dthalf = 0.5 * update->dt;
dt4 = 0.25 * update->dt;
dt8 = 0.125 * update->dt;
tdrag_factor = 1.0 - (update->dt * t_freq * drag);
}
double FixNVTManifoldRattle::memory_usage()
{
double bytes = FixNVEManifoldRattle::memory_usage();
bytes += (4*mtchain+1)*sizeof(double);
return bytes;
}

View File

@ -0,0 +1,152 @@
/* ----------------------------------------------------------------------
Lammps - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
-----------------------------------------------------------------------
This file is a part of the USER-MANIFOLD package.
Copyright (2013-2014) Stefan Paquay, Eindhoven University of Technology.
License: GNU General Public License.
See the README file in the top-level LAMMPS directory.
This file is part of the user-manifold package written by
Stefan Paquay at the Eindhoven University of Technology.
This module makes it possible to do MD with particles constrained
to pretty arbitrary manifolds characterised by some constraint function
g(x,y,z) = 0 and its normal grad(g). The number of manifolds available
right now is limited but can be extended straightforwardly by making
a new class that inherits from manifold and implements all pure virtual
methods.
Thanks to Remy Kusters for beta-testing!
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(nvt/manifold/rattle,FixNVTManifoldRattle)
#else
#ifndef LMP_FIX_NVT_MANIFOLD_RATTLE_H
#define LMP_FIX_NVT_MANIFOLD_RATTLE_H
#include "fix_nve_manifold_rattle.h"
/*
FixNVTManifoldRattle works by wrapping some Nose-Hoover thermostat routines
around the time integration functions of FixNVEManifoldRattle.
*/
namespace LAMMPS_NS {
class FixNVTManifoldRattle : public FixNVEManifoldRattle {
public:
FixNVTManifoldRattle(LAMMPS *, int, char **, int = 1);
virtual ~FixNVTManifoldRattle();
virtual void initial_integrate(int);
virtual void final_integrate();
virtual void init();
virtual void reset_dt();
virtual int setmask();
virtual void setup(int); // Not needed for fixNVE but is for fixNVT
virtual double memory_usage();
protected:
void compute_temp_target();
void nhc_temp_integrate();
void nh_v_temp();
double dthalf, dt4, dt8;
char *id_temp;
class Compute* temperature;
double t_start,t_stop, t_period;
double t_current,t_target,ke_target;
double t_freq, drag, tdrag_factor;
double boltz,nktv2p,tdof;
double *eta,*eta_dot;
double *eta_dotdot;
double *eta_mass;
int mtchain;
double factor_eta;
int which, got_temp;
const char *fix_id;
};
}
#endif // LMP_FIX_NVE_MANIFOLD_RATTLE_H
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: There is no manifold named ...
Self-explanatory. You requested a manifold whose name was not
registered at the factory.
E: Manifold pointer was NULL for some reason!
This indicates a bug. The factory was unable to properly create
the requested manifold even though it was registered. Send the
maintainer an e-mail.
E: Manifold ... needs at least ... argument(s)!
Self-explanatory. Provide enough arguments for the proper
creating of the requested manifold.
E: Parameter pointer was NULL!
This indicates a bug. The array that contains the parameters
could not be allocated. Send the maintainer an e-mail.
E: Could not allocate space for arg!
One of the arguments provided was too long (it contained
too many characters)
E: Option ... needs ... argument(s):
Self-explanatory. Read the documentation of this fix properly.
E: Illegal fix nve/manifold/rattle command! Option ... not recognized!
Self-explanatory. The given option(s) do not exist.
E: Variable name for fix nve/manifold/rattle does not exist
Self-explanatory.
E: Variable for fix nve/manifold/rattle is invalid style
fix nve/manifold/rattle only supports equal style variables.
*/

View File

@ -0,0 +1,101 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
-----------------------------------------------------------------------
This file is a part of the USER-MANIFOLD package.
This package allows LAMMPS to perform MD simulations of particles
constrained on a manifold (i.e., a 2D subspace of the 3D simulation
box). It achieves this using the RATTLE constraint algorithm applied
to single-particle constraint functions g(xi,yi,zi) = 0 and their
derivative (i.e. the normal of the manifold) n = grad(g).
It is very easy to add your own manifolds to the current zoo
(we now have sphere, a dendritic spine approximation, a 2D plane (for
testing purposes) and a wave-y plane.
See the README file for more info.
Stefan Paquay, s.paquay@tue.nl
Applied Physics/Theory of Polymers and Soft Matter,
Eindhoven University of Technology (TU/e), The Netherlands
Thanks to Remy Kusters at TU/e for testing.
This software is distributed under the GNU General Public License.
------------------------------------------------------------------------- */
#ifndef LMP_MANIFOLD_H
#define LMP_MANIFOLD_H
#include "pointers.h"
#include "error.h"
#include "math.h"
#include "string.h"
#include "lmptype.h"
namespace LAMMPS_NS {
// Abstract base class.
class manifold {
public:
manifold() : params(NULL){ }
virtual ~manifold(){ delete[] params; }
virtual double g( const double * ) = 0;
virtual void n( const double *, double * ) = 0;
// Variant of g that computes n at the same time.
virtual double g_and_n( const double *x, double *nn )
{
this->n(x,nn);
return g(x);
}
virtual const char *id() = 0;
virtual void set_atom_id( tagint a_id ){}
virtual int nparams() = 0;
double **get_params(){ return &params; };
// Overload if any initialisation depends on params:
virtual void post_param_init(){}
virtual void checkup(){} // Some diagnostics...
protected:
double *params;
};
// Some utility functions that are templated, so I implement them
// here in the header.
template< unsigned int size > inline
double infnorm( double *vect )
{
double largest = fabs( vect[0] );
for( unsigned int i = 1; i < size; ++i ){
double c = fabs( vect[i] );
largest = ( c > largest ) ? c : largest;
}
return largest;
}
inline double dot( double *a, double *b ){
return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
}
} // namespace LAMMPS_NS
#endif // LMP_MANIFOLD_H

View File

@ -0,0 +1,22 @@
#include "manifold_cylinder.h"
using namespace LAMMPS_NS;
double manifold_cylinder::g( const double *x )
{
double R = params[0];
double r2 = x[0]*x[0] + x[1]*x[1];
return R*R - r2;
}
void manifold_cylinder::n( const double *x, double *n )
{
n[0] = -2*x[0];
n[1] = -2*x[1];
n[2] = 0.0;
}
manifold_cylinder::manifold_cylinder( LAMMPS *lmp, int argc,
char **argv ) : Pointers(lmp)
{}

View File

@ -0,0 +1,26 @@
#ifndef LMP_MANIFOLD_CYLINDER_H
#define LMP_MANIFOLD_CYLINDER_H
#include "manifold.h"
// A normal cylinder
namespace LAMMPS_NS {
class manifold_cylinder : public manifold, protected Pointers {
public:
enum { NPARAMS = 1 }; // Number of parameters.
manifold_cylinder( LAMMPS *lmp, int, char ** );
virtual ~manifold_cylinder(){}
virtual double g( const double *x );
virtual void n( const double *x, double *n );
static const char *type(){ return "cylinder"; }
virtual const char *id(){ return type(); }
static int expected_argc(){ return NPARAMS; }
virtual int nparams(){ return NPARAMS; }
};
}
#endif // LMP_MANIFOLD_CYLINDER_H

View File

@ -0,0 +1,42 @@
#include "manifold_cylinder_dent.h"
#include "math_const.h"
#include <math.h>
using namespace LAMMPS_NS;
manifold_cylinder_dent::manifold_cylinder_dent( LAMMPS *lmp, int argc,
char **argv ) : Pointers(lmp)
{}
double manifold_cylinder_dent::g( const double *x )
{
double l = params[1], R = params[0], a = params[2];
double r2 = x[1]*x[1] + x[0]*x[0];
if( fabs(x[2]) < 0.5*l ){
double k = MathConst::MY_2PI / l;
double c = R - 0.5*a*( 1.0 + cos(k*x[2]) );
return c*c - r2;
}else{
return R*R - r2;
}
}
void manifold_cylinder_dent::n( const double *x, double *n )
{
double l = params[1], R = params[0], a = params[2];
if( fabs(x[2]) < 0.5*l ){
double k = MathConst::MY_2PI / l;
double c = R - 0.5*a*(1.0 + cos(k*x[2]));
n[0] = -2*x[0];
n[1] = -2*x[1];
n[2] = c*a*k*sin(k*x[2]);
}else{
n[0] = -2*x[0];
n[1] = -2*x[1];
n[2] = 0.0;
}
}

View File

@ -0,0 +1,24 @@
#ifndef LMP_MANIFOLD_CYLINDER_DENT_H
#define LMP_MANIFOLD_CYLINDER_DENT_H
#include "manifold.h"
namespace LAMMPS_NS {
class manifold_cylinder_dent : public manifold, protected Pointers {
public:
manifold_cylinder_dent( LAMMPS *lmp, int, char ** );
enum { NPARAMS = 3 }; // Number of parameters.
virtual ~manifold_cylinder_dent(){}
virtual double g( const double *x );
virtual void n( const double *x, double *n );
static const char *type(){ return "cylinder/dent"; }
virtual const char *id(){ return type(); }
static int expected_argc(){ return NPARAMS; }
virtual int nparams(){ return NPARAMS; }
};
}
#endif // LMP_MANIFOLD_CYLINDER_DENT_H

View File

@ -0,0 +1,41 @@
#include "manifold_dumbbell.h"
#include <math.h>
using namespace LAMMPS_NS;
manifold_dumbbell::manifold_dumbbell( LAMMPS *lmp, int argc, char **argv ) : Pointers(lmp)
{}
/*
* Equation for dumbbell is:
*
* -x^2 - y^2 + (a^2 - (z/c)^2)*( 1 + A*sin(const(z) *z^2) )^4
* const(z) = (z < 0) ? B2 : B
* params[4] = { a, A, B, c }
*/
double manifold_dumbbell::g( const double *x )
{
double a = params[0];
double A = params[1];
double B = params[2];
double c = params[3];
return -1.0*(x[0]*x[0]+x[1]*x[1])+(a*a-x[2]*x[2]/(c*c))*(1.0+pow(A*sin(B *x[2]*x[2]),4));
}
void manifold_dumbbell::n( const double *x, double *nn )
{
double a = params[0];
double A = params[1];
double B = params[2];
double c = params[3];
nn[0] = -2.0*x[0];
nn[1] = -2.0*x[1];
nn[2] = 8.0*A*A*A*A*B*x[2]*(a*a-(x[2]*x[2]/(c*c)))*cos(B*x[2]*x[2]) *
pow(sin(B*x[2]*x[2]),3)-2*x[2]*(1+pow(A*sin(B*x[2]*x[2]),4))/(c*c);
}

View File

@ -0,0 +1,28 @@
#ifndef LMP_MANIFOLD_DUMBBELL_H
#define LMP_MANIFOLD_DUMBBELL_H
#include "manifold.h"
namespace LAMMPS_NS {
// A dendritic dumbbell approximation:
class manifold_dumbbell : public manifold, protected Pointers {
public:
enum { NPARAMS = 4 }; // Number of parameters.
manifold_dumbbell( LAMMPS *lmp, int, char ** );
virtual ~manifold_dumbbell(){}
virtual double g ( const double *x );
virtual void n ( const double *x, double *nn );
static const char* type(){ return "dumbbell"; }
virtual const char *id(){ return type(); }
static int expected_argc(){ return NPARAMS; }
virtual int nparams(){ return NPARAMS; }
};
}
#endif // LMP_MANIFOLD_DUMBBELL_H

View File

@ -0,0 +1,26 @@
#include "manifold_ellipsoid.h"
using namespace LAMMPS_NS;
manifold_ellipsoid::manifold_ellipsoid( LAMMPS *lmp, int narg, char **argv ) : Pointers(lmp)
{}
double manifold_ellipsoid::g( const double *x )
{
const double ai = 1.0 / params[0];
const double bi = 1.0 / params[1];
const double ci = 1.0 / params[2];
return x[0]*x[0]*ai*ai + x[1]*x[1]*bi*bi + x[2]*x[2]*ci*ci - 1.0;
}
void manifold_ellipsoid::n( const double *x, double * n )
{
const double ai = 1.0 / params[0];
const double bi = 1.0 / params[1];
const double ci = 1.0 / params[2];
n[0] = 2*x[0]*ai*ai;
n[1] = 2*x[1]*bi*bi;
n[2] = 2*x[2]*ci*ci;
}

View File

@ -0,0 +1,26 @@
#ifndef LMP_MANIFOLD_ELLIPSOID_H
#define LMP_MANIFOLD_ELLIPSOID_H
#include "manifold.h"
namespace LAMMPS_NS {
// An ellipsoid:
class manifold_ellipsoid : public manifold, protected Pointers {
public:
enum { NPARAMS = 3 };
manifold_ellipsoid( LAMMPS *lmp, int, char ** );
virtual ~manifold_ellipsoid(){}
virtual double g( const double *x );
virtual void n( const double *x, double *n );
static const char* type(){ return "ellipsoid"; }
virtual const char *id(){ return type(); }
static int expected_argc(){ return NPARAMS; }
virtual int nparams(){ return NPARAMS; }
};
}
#endif // LMP_MANIFOLD_ELLIPSOID_H

View File

@ -0,0 +1,154 @@
/* ----------------------------------------------------------------------
Lammps - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
-----------------------------------------------------------------------
This file is a part of the USER-MANIFOLD package.
Copyright (2013-2014) Stefan Paquay, Eindhoven University of Technology.
License: GNU General Public License.
See the README file in the top-level LAMMPS directory.
This file is part of the user-manifold package written by
Stefan Paquay at the Eindhoven University of Technology.
This module makes it possible to do MD with particles constrained
to pretty arbitrary manifolds characterised by some constraint function
g(x,y,z) = 0 and its normal grad(g). The number of manifolds available
right now is limited but can be extended straightforwardly by making
a new class that inherits from manifold and implements all pure virtual
methods.
Thanks to Remy Kusters for beta-testing!
------------------------------------------------------------------------- */
#ifndef LMP_MANIFOLD_FACTORY_H
#define LMP_MANIFOLD_FACTORY_H
#include "manifold.h"
#include "pointers.h"
#include <string>
#include <map>
#include "manifold_cylinder.h"
#include "manifold_cylinder_dent.h"
#include "manifold_dumbbell.h"
#include "manifold_ellipsoid.h"
#include "manifold_plane.h"
#include "manifold_plane_wiggle.h"
#include "manifold_sphere.h"
#include "manifold_supersphere.h"
#include "manifold_spine.h"
#include "manifold_spine_opt.h"
#include "manifold_thylakoid.h"
#include "manifold_torus.h"
inline int str_eql(const char *s1, const char *s2)
{
return strcmp(s1,s2) == 0;
}
/*
* Defining USE_PHONY_LAMMPS makes sure that none of the LAMMPS classes are
* included/compiled. This is done in order to allow other programs to use
* the manifold_factory without compiling all of LAMMPS itself. The relevant
* classes/functions are replaced with dummy ones defined in this #ifdef-block:
*/
#ifdef USE_PHONY_LAMMPS
# ifdef __GNUG__
# warning Not compiling actual LAMMPS classes!
# endif
struct Error {
void all(const char *fname, int line, const char* msg)
{
fprintf(stderr,"ERROR: %s (%s:%d)",msg,fname,line);
std::terminate();
}
void one(const char *fname, int line, const char* msg)
{
fprintf(stderr,"ERROR: %s (%s:%d)",msg,fname,line);
std::terminate();
}
};
struct LAMMPS { };
struct Pointers
{
Pointers(LAMMPS *) : error( &e ){}
Error e;
Error *error;
};
static FILE *screen = fopen("/dev/stdout","w");
#define FLERR __FILE__,__LINE__ // Equivalent to definition in pointers.h
#endif // USE_PHONY_LAMMPS
// Macro to simplify the fix_impl creation:
#define RETURN_MANIFOLD_IF(MNAME, MANIFOLD_TYPE) \
do { \
using namespace LAMMPS_NS; \
if( strcmp( MNAME, MANIFOLD_TYPE::type() ) == 0 ) {\
return new MANIFOLD_TYPE(lmp,narg,arg); \
} \
}while(0)
namespace LAMMPS_NS {
template <typename m_type>
void make_manifold_if( manifold **man_ptr, const char *name,
LAMMPS *lmp, int narg, char **arg )
{
if( strcmp( m_type::type(), name ) == 0 ){
if( *man_ptr == NULL ){
*man_ptr = new m_type(lmp, narg, arg);
}
}
}
inline manifold* create_manifold(const char *mname, LAMMPS *lmp,
int narg, char **arg )
{
using namespace LAMMPS_NS;
manifold *man = NULL;
make_manifold_if<manifold_cylinder> ( &man, mname, lmp, narg, arg );
make_manifold_if<manifold_cylinder_dent>( &man, mname, lmp, narg, arg );
make_manifold_if<manifold_dumbbell> ( &man, mname, lmp, narg, arg );
make_manifold_if<manifold_ellipsoid> ( &man, mname, lmp, narg, arg );
make_manifold_if<manifold_plane> ( &man, mname, lmp, narg, arg );
make_manifold_if<manifold_plane_wiggle> ( &man, mname, lmp, narg, arg );
make_manifold_if<manifold_sphere> ( &man, mname, lmp, narg, arg );
make_manifold_if<manifold_supersphere> ( &man, mname, lmp, narg, arg );
make_manifold_if<manifold_spine> ( &man, mname, lmp, narg, arg );
// make_manifold_if<manifold_spine_opt> ( &man, mname, lmp, narg, arg );
make_manifold_if<manifold_thylakoid> ( &man, mname, lmp, narg, arg );
make_manifold_if<manifold_torus> ( &man, mname, lmp, narg, arg );
return man;
}
} // namespace LAMMPS_NS
#endif // LMP_MANIFOLD_FACTORY_H

View File

@ -0,0 +1,54 @@
#ifndef LMP_MANIFOLD_SUPERSPHERE_H
#define LMP_MANIFOLD_SUPERSPHERE_H
#include "manifold.h"
namespace LAMMPS_NS {
// A supersphere:
class manifold_supersphere : public manifold, protected Pointers {
public:
enum { NPARAMS = 2 };
manifold_supersphere( LAMMPS *lmp, int, char ** ) : Pointers(lmp){}
virtual ~manifold_supersphere(){}
double my_sign( double a )
{
return (a > 0) - (a < 0);
}
virtual double g( const double *x )
{
double R = params[0];
double q = params[1];
double xx = fabs(x[0]);
double yy = fabs(x[1]);
double zz = fabs(x[2]);
double rr = pow(xx,q) + pow(yy,q) + pow(zz,q);
return rr - pow(R,q);
}
virtual void n( const double *x, double *nn )
{
double q = params[1];
double xx = fabs(x[0]);
double yy = fabs(x[1]);
double zz = fabs(x[2]);
nn[0] = q * my_sign(x[0])*pow(xx,q-1);
nn[1] = q * my_sign(x[1])*pow(yy,q-1);
nn[2] = q * my_sign(x[2])*pow(zz,q-1);
}
static const char* type(){ return "supersphere"; }
virtual const char *id(){ return type(); }
static int expected_argc(){ return NPARAMS; }
virtual int nparams(){ return NPARAMS; }
};
}
#endif

View File

@ -0,0 +1,22 @@
#include "manifold_plane.h"
using namespace LAMMPS_NS;
manifold_plane::manifold_plane( LAMMPS *lmp, int argc, char **argv ) :
Pointers(lmp)
{}
double manifold_plane::g( const double *x )
{
double a = params[0], b = params[1], c = params[2];
double x0 = params[3], y0 = params[4], z0 = params[5];
return a*(x[0] - x0) + b*(x[1] - y0) + c*(x[2] - z0);
}
void manifold_plane::n( const double *x, double *n )
{
n[0] = params[0];
n[1] = params[1];
n[2] = params[2];
}

View File

@ -0,0 +1,25 @@
#ifndef LMP_MANIFOLD_PLANE_H
#define LMP_MANIFOLD_PLANE_H
#include "manifold.h"
namespace LAMMPS_NS {
// A 2D plane
class manifold_plane : public manifold, protected Pointers {
public:
enum { NPARAMS = 6 }; // Number of parameters.
manifold_plane( LAMMPS *lmp, int, char ** );
virtual ~manifold_plane(){}
virtual double g( const double *x );
virtual void n( const double *x, double *n );
static const char *type(){ return "plane"; }
virtual const char *id(){ return type(); }
static int expected_argc(){ return NPARAMS; }
virtual int nparams(){ return NPARAMS; }
};
}
#endif // LMP_MANIFOLD_PLANE_H

View File

@ -0,0 +1,27 @@
#include "manifold_plane_wiggle.h"
#include <math.h>
using namespace LAMMPS_NS;
manifold_plane_wiggle::manifold_plane_wiggle( LAMMPS *lmp, int argc, char **argv ) :
Pointers(lmp)
{}
double manifold_plane_wiggle::g( const double *x )
{
double a = params[0];
double w = params[1];
return x[2] - a*sin(w*x[0]);
}
void manifold_plane_wiggle::n( const double *x, double *n )
{
double a = params[0];
double w = params[1];
n[2] = 1;
n[1] = 0.0;
n[0] = -a*w*cos(x[0]);
}

View File

@ -0,0 +1,25 @@
#ifndef LMP_MANIFOLD_PLANE_WIGGLE_H
#define LMP_MANIFOLD_PLANE_WIGGLE_H
#include "manifold.h"
namespace LAMMPS_NS {
// A 2D wiggly/wave-y plane (Like z = A cos(kx))
class manifold_plane_wiggle : public manifold, protected Pointers {
public:
enum { NPARAMS = 2 }; // Number of parameters.
manifold_plane_wiggle( LAMMPS *lmp, int, char ** );
virtual ~manifold_plane_wiggle(){}
virtual double g( const double *x );
virtual void n( const double *x, double *n );
static const char *type(){ return "plane/wiggle"; }
virtual const char *id(){ return type(); }
static int expected_argc(){ return NPARAMS; }
virtual int nparams(){ return NPARAMS; }
};
}
#endif // LMP_MANIFOLD_PLANE_WIGGLE_H

View File

@ -0,0 +1,695 @@
/* ----------------------------------------------------------------------
Lammps - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
-----------------------------------------------------------------------
This file is a part of the USER-MANIFOLD package.
Copyright (2013-2014) Stefan Paquay, Eindhoven University of Technology.
License: GNU General Public License.
See the README file in the top-level LAMMPS directory.
This file is part of the user-manifold package written by
Stefan Paquay at the Eindhoven University of Technology.
This module makes it possible to do MD with particles constrained
to pretty arbitrary manifolds characterised by some constraint function
g(x,y,z) = 0 and its normal grad(g). The number of manifolds available
right now is limited but can be extended straightforwardly by making
a new class that inherits from manifold and implements all pure virtual
methods.
Thanks to Remy Kusters for beta-testing!
------------------------------------------------------------------------- */
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "input.h"
#include "manifold_factory.h"
#include "update.h"
#include "variable.h"
#include "manifold_rattle_integrator.h"
using namespace LAMMPS_NS;
enum { CONST, EQUAL }; // For treating the variables.
inline int was_var( const char *str )
{
if( strlen(str) > 2 ){
return (str[0] == 'v') && (str[1] == '_');
}else{
return 0;
}
}
void LAMMPS_NS::shift_args_back( int argi, char ***arg_ptr, int delta, int &narg )
{
for( int i = argi; i < narg-delta; ++i ){
(*arg_ptr)[i] = (*arg_ptr[i+delta]);
}
narg -= delta;
}
/* -------------------------------------------------------------------------
Integrator stuff. Constructor:
------------------------------------------------------------------------ */
ManifoldRattleIntegratorBase::ManifoldRattleIntegratorBase(LAMMPS *lmp, int &narg,
char ***arg_ptr )
: Pointers(lmp), tstrs(NULL), tvars(NULL), tstyle(NULL), is_var(NULL)
{
char **arg = *arg_ptr;
// Defaults for member variables:
nevery = 0;
dtv = dtf = 0;
tolerance = force->numeric( FLERR, arg[3] );
max_iter = force->numeric( FLERR, arg[4] );
update_style = 1;
int argi = 6; // 6 for general args, + nvars for manifold vars.
while( argi < narg ){
if( comm->me == 0 ){
fprintf(screen,"Parsing arg \"%s\"\n",arg[argi]);
}
if( is_legal_keyword( arg[argi] ) ){
int status = parse_arg(argi, narg, arg_ptr);
if( status ){
char msg[256];
sprintf(msg,"Error parsing arg \"%s\".\n", arg[argi]);
error->all(FLERR, msg);
}
}else{
++argi;
}
}
}
/* -------------------------------------------------------------------------
Destructor:
------------------------------------------------------------------------ */
ManifoldRattleIntegratorBase::~ManifoldRattleIntegratorBase()
{
if( tstrs ){
for( int i = 0; i < nvars; ++i ){
delete [] tstrs[i];
}
delete [] tstrs;
}
if( tvars ) delete [] tvars;
if( tstyle ) delete [] tstyle;
if( is_var ) delete [] is_var;
}
void ManifoldRattleIntegratorBase::init_manifold( int narg, char **arg,
manifold *man_ptr )
{
if( !man_ptr ) error->all(FLERR,"manifold ptr was NULL");
tstrs = new char*[nvars];
tvars = new int[nvars];
tstyle = new int[nvars];
is_var = new int[nvars];
if( !tstrs || !tvars || !tstyle || !is_var ){
error->all(FLERR, "Error creating manifold arg arrays");
}
for( int i = 0; i < nvars; ++i ){
int len = 0, offset = 0;
if( was_var( arg[i+6] ) ){
len = strlen(arg[i+6]) - 1; // -1 because -2 for v_, +1 for \0.
is_var[i] = 1;
offset = 2;
}else{
force->numeric(FLERR,arg[i+6]); // Check if legal number.
len = strlen( arg[i+6] ) + 1; // +1 for \0.
is_var[i] = 0;
}
tstrs[i] = new char[len];
if( tstrs[i] == NULL ) error->all(FLERR,"Error allocating space for args.");
strcpy( tstrs[i], arg[i+6] + offset );
}
*(man_ptr->get_params()) = new double[nvars];
if( !(*man_ptr->get_params()) ) error->all(FLERR,"Failed to allocate params!");
for( int i = 0; i < nvars; ++i ){
// If param i was variable type, it will be set later...
(*man_ptr->get_params())[i] = is_var[i] ? 0.0 : force->numeric( FLERR, arg[i+6] );
}
man_ptr->post_param_init();
}
void ManifoldRattleIntegratorBase::reset_dt()
{
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
}
// Identifies keyword in (*arg_ptr)[argi], handles it properly, and strips
// the relevant args from (*arg_ptr) and subtracts the number of args handled
// from narg (so that (*arg_ptr) contains narg elements again.
// Returns 0 on success, -1 otherwise.
int ManifoldRattleIntegratorBase::parse_arg( int argi, int &narg,
char ***arg_ptr )
{
if( strcmp((*arg_ptr)[argi], "every") == 0 ){
if( argi + 1 >= narg ){
error->all(FLERR,"keyword every needs one argument");
}
nevery = force->inumeric(FLERR,(*arg_ptr)[argi+1]);
shift_args_back( argi, arg_ptr, 2, narg );
return 0;
}
if( strcmp((*arg_ptr)[argi], "style") == 0 ){
if( argi + 1 >= narg ){
error->all(FLERR,"keyword every needs one argument");
}
update_style = force->inumeric(FLERR, (*arg_ptr)[argi+1]);
if( update_style < 1 || update_style > 2 ){
error->warning(FLERR,"update style not known, reverting to standard");
}
shift_args_back( argi, arg_ptr, 2, narg );
return 0;
}
return -1;
}
int ManifoldRattleIntegratorBase::is_legal_keyword( const char *arg )
{
return strcmp(arg, "every") == 0;
}
void ManifoldRattleIntegratorBase::print_stats(const char *header)
{
double n = stats.natoms;
if( n > 0 ){
stats.x_iters_per_atom += stats.x_iters / n;
stats.v_iters_per_atom += stats.v_iters / n;
}
double x_iters = 0, v_iters = 0;
bigint ntimestep = update->ntimestep;
int me = -1;
MPI_Comm_rank(world,&me);
MPI_Allreduce(&stats.x_iters_per_atom,&x_iters,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&stats.v_iters_per_atom,&v_iters,1,MPI_DOUBLE,MPI_SUM,world);
// Set iters back to zero:
stats.x_iters_per_atom = stats.x_iters = 0;
stats.v_iters_per_atom = stats.v_iters = 0;
if( me == 0 ){
double inv_tdiff = 1.0/( static_cast<double>(ntimestep) - stats.last_out );
stats.last_out = ntimestep;
fprintf(screen, "%s stats for time step " BIGINT_FORMAT " on %d atoms:\n",
header, ntimestep, stats.natoms);
fprintf(screen, " iters/atom: x = %f, v = %f, dofs removed %d",
x_iters * inv_tdiff, v_iters * inv_tdiff, stats.dofs_removed);
fprintf(screen,"\n");
}
}
/* ----------------------------------------------------------------------------
Modify the force so that the new positions remain approximately on manifold.
Time step velocity (half-way) and position accordingly.
-------------------------------------------------------------------------------- */
void ManifoldRattleIntegratorBase::nve_x_rattle(int igroup, int groupbit)
{
double dtfm;
// update v and x of atoms in group
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int natoms = 0;
if (igroup == atom->firstgroup){
nlocal = atom->nfirst;
}
if (rmass) {
for (int i = 0; i < nlocal; i++){
if (mask[i] & groupbit){
natoms++;
dtfm = dtf / rmass[i];
rattle_manifold_x( x[i], v[i], f[i], dtv, dtfm, atom->tag[i] );
}
}
} else {
for (int i = 0; i < nlocal; i++){
if (mask[i] & groupbit) {
natoms++;
dtfm = dtf / mass[type[i]];
rattle_manifold_x( x[i], v[i], f[i], dtv, dtfm, atom->tag[i] );
}
}
}
if( nevery > 0 ){
// Count ALL atoms this fix works on:
MPI_Allreduce(&natoms,&stats.natoms,1,MPI_INT,MPI_SUM,world);
}
}
/* ----------------------------------------------------------------------------
Time step (other half) the velocities so that they remain in manifold.
-------------------------------------------------------------------------------- */
void ManifoldRattleIntegratorBase::nve_v_rattle(int igroup, int groupbit)
{
double dtfm;
// update v of atoms in group
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
rattle_manifold_v( v[i], f[i], x[i], dtfm );
}
}
} else {
for (int i = 0; i < nlocal; i++){
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
rattle_manifold_v( v[i], f[i], x[i], dtfm );
}
}
}
}
/******************* ManifoldRattleIntegrator stuff. *****************/
/* ----------------------------------------------------------------------
Performs a rattle update for x
------------------------------------------------------------------------- */
void ManifoldRattleIntegrator::rattle_manifold_x(double *x, double *v,
double *f, double dtv,
double dtfm, int tagi )
{
/*
A RATTLE update for the position constraint.
Original update is x += dtv * v_1/2
Now you do
v_1/2(lambda) = v_0 + dtfm * ( f + lambda*n_old )
and solve
xold - xnew + dtv * v_1/2(lambda) = 0
g(xnew) = 0
for x and lambda. The lambda you find then gives v_1/2 as well.
*/
double xo[3]; // Previous position to update from.
double vo[3]; // Previous velocity to update from.
double l = 0; // Lagrangian multiplier for constraint forces.
double R[4]; // System that is 0.
double dx[4]; // Update that follows from Newton iteration.
double no[3]; // Normal at xo.
double nn[3]; // Normal at x, the new position.
double res; // Residual.
int iters = 0; // Iterations used
double c = dtfm*dtv; // Used for iterating in the Newton loop:
double no_nn, nn_R;
vo[0] = v[0];
vo[1] = v[1];
vo[2] = v[2];
xo[0] = x[0];
xo[1] = x[1];
xo[2] = x[2];
double gg = ptr_m->g_and_n(x,no);
nn[0] = no[0];
nn[1] = no[1];
nn[2] = no[2];
double vt[3];
vt[0] = vo[0] + dtfm*f[0];
vt[1] = vo[1] + dtfm*f[1];
vt[2] = vo[2] + dtfm*f[2];
double no_dt[3];
no_dt[0] = dtfm*no[0];
no_dt[1] = dtfm*no[1];
no_dt[2] = dtfm*no[2];
// Assume that no_nn is roughly constant during iteration:
const double c_inv = 1.0 / c;
while ( 1 ) {
v[0] = vt[0] - l*no_dt[0];
v[1] = vt[1] - l*no_dt[1];
v[2] = vt[2] - l*no_dt[2];
R[0] = xo[0] - x[0] + dtv * v[0];
R[1] = xo[1] - x[1] + dtv * v[1];
R[2] = xo[2] - x[2] + dtv * v[2];
R[3] = gg;
// Analytic solution to system J*(dx,dy,dz,dl)^T = R
// no_nn = no[0]*nn[0] + no[1]*nn[1] + no[2]*nn[2];
nn_R = nn[0]*R[0] + nn[1]*R[1] + nn[2]*R[2];
no_nn = no[0]*nn[0] + no[1]*nn[1] + no[2]*nn[2];
double n_inv = 1.0 / no_nn;
// fprintf( screen, "nn_R = %f, no_nn = %f\n", nn_R, no_nn );
dx[3] = -nn_R - R[3];
dx[3] *= n_inv;
dx[0] = -R[0] - no[0]*dx[3];
dx[1] = -R[1] - no[1]*dx[3];
dx[2] = -R[2] - no[2]*dx[3];
dx[3] *= c_inv;
double smooth = 1.0; // / sqrt( 1.0 + res*res );
x[0] -= smooth*dx[0];
x[1] -= smooth*dx[1];
x[2] -= smooth*dx[2];
l -= smooth*dx[3];
/*
fprintf( screen, "dx = (%e, %e, %e, %e), R = (%e, %e, %e, %e)\n",
dx[0], dx[1], dx[2], dx[3], R[0], R[1], R[2], R[3] );
fprintf( screen, "x = (%f, %f, %f, %f), n = (%f, %f, %f)\n",
x[0], x[1], x[2], l, nn[0], nn[1], nn[2] );
*/
res = infnorm<4>(R);
++iters;
if( (res < tolerance) || (iters >= max_iter) ) break;
// Update nn and g.
gg = ptr_m->g(x);
ptr_m->n(x,nn);
// gg = ptr_m->g(x);
}
if( iters >= max_iter && res > tolerance ){
char msg[2048];
sprintf(msg,"Failed to constrain atom %d (x = (%f, %f, %f)! res = %e, iters = %d\n",
tagi, x[0], x[1], x[2], res, iters);
error->one(FLERR,msg);
}
// "sync" x and v:
v[0] = vt[0] - l*no_dt[0];
v[1] = vt[1] - l*no_dt[1];
v[2] = vt[2] - l*no_dt[2];
stats.x_iters += iters;
}
/* ----------------------------------------------------------------------
Performs a rattle update for v
------------------------------------------------------------------------- */
void ManifoldRattleIntegrator::rattle_manifold_v(double *v, double *f,
double *x, double dtfm,
int tagi )
{
/*
The original update was
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
Now you add the rattle-like update:
vold - vnew + dtfm * F + mu * n_new = 0
dot( vnew, n_new ) = 0
*/
double vo[3]; // V at t + 1/2 dt
double l = 0; // Lagrangian multiplier for constraint forces.
double R[4]; // System that is 0.
double dv[4]; // Update that follows from Newton iteration.
double n[3]; // Normal.
double res; // Residual.
int iters = 0; // Iterations used
double c = dtfm; // Used for iterating in the Newton loop:
double nn2, nn_R;
vo[0] = v[0];
vo[1] = v[1];
vo[2] = v[2];
// Initial guess is unconstrained update:
v[0] += dtfm*f[0];
v[1] += dtfm*f[1];
v[2] += dtfm*f[2];
ptr_m->n(x,n);
double vt[3];
vt[0] = vo[0] + dtfm*f[0];
vt[1] = vo[1] + dtfm*f[1];
vt[2] = vo[2] + dtfm*f[2];
double no_dt[3];
no_dt[0] = dtfm*n[0];
no_dt[1] = dtfm*n[1];
no_dt[2] = dtfm*n[2];
nn2 = n[0]*n[0] + n[1]*n[1] + n[2]*n[2];
const double n_inv = 1.0 / nn2;
const double c_inv = 1.0 / c;
switch( update_style ){
default:
case 1:
do{
R[0] = vt[0] - v[0] - l * no_dt[0];
R[1] = vt[1] - v[1] - l * no_dt[1];
R[2] = vt[2] - v[2] - l * no_dt[2];
R[3] = v[0]*n[0] + v[1]*n[1] + v[2]*n[2];
// Analytic solution to system J*(dx,dy,dz,dl)^T = R
nn_R = n[0]*R[0] + n[1]*R[1] + n[2]*R[2];
dv[3] = -nn_R - R[3];
dv[3] *= n_inv;
dv[0] = -n[0]*dv[3] - R[0];
dv[1] = -n[1]*dv[3] - R[1];
dv[2] = -n[2]*dv[3] - R[2];
dv[3] *= c_inv;
v[0] -= dv[0];
v[1] -= dv[1];
v[2] -= dv[2];
l -= dv[3];
res = infnorm<4>(R);
++iters;
}while( (res > tolerance) && (iters < max_iter) );
if( iters >= max_iter && res >= tolerance ){
char msg[2048];
sprintf(msg,"Failed to constrain atom %d (x = (%f, %f, %f)! res = %e, iters = %d\n",
tagi, x[0], x[1], x[2], res, iters);
error->all(FLERR,msg);
}
break;
case 2:{
double mass = dtfm / update->dt;
double vhnp = vo[0]*n[0] + vo[1]*n[1] + vo[2]*n[2];
double fpnp = f[0]*n[0] + f[1]*n[1] + f[2]*n[2];
double phnp = dtfm*vhnp / update->dt;
l = (2.0*mass/update->dt)*phnp + fpnp;
v[0] += 0.5*dtv*f[0] - n[0]*(phnp + fpnp);
v[1] += 0.5*dtv*f[1] - n[1]*(phnp + fpnp);
v[2] += 0.5*dtv*f[2] - n[2]*(phnp + fpnp);
break;
}
}
stats.v_iters += iters;
}
int ManifoldRattleIntegratorBase::dof( int igroup, int groupbit )
{
int *mask = atom->mask;
int nlocal = atom->nlocal;
int natoms = 0;
for( int i = 0; i < nlocal; ++i ){
if(mask[i] & groupbit) ++natoms;
}
int dofs;
MPI_Allreduce( &natoms, &dofs, 1, MPI_INT, MPI_SUM, world );
// Make sure that, if there is just no or one atom, no dofs are subtracted,
// since for the first atom already 3 dofs are subtracted because of the
// centre of mass corrections:
if( dofs <= 1 ) dofs = 0;
stats.dofs_removed = dofs;
return dofs;
}
ManifoldRattleIntegrator::ManifoldRattleIntegrator( LAMMPS *lmp, int &narg,
char ***arg_ptr ) :
ManifoldRattleIntegratorBase( lmp, narg, arg_ptr ), ptr_m(NULL)
{
char **arg = *arg_ptr;
ptr_m = create_manifold(arg[5], lmp, narg, arg);
if( !ptr_m ){
char msg[2048];
sprintf(msg,"Could not create manifold of type \"%s\"", arg[5]);
error->all(FLERR,msg);
}
nvars = ptr_m->nparams();
init_manifold( narg, arg, ptr_m );
}
ManifoldRattleIntegrator::~ManifoldRattleIntegrator()
{
if( ptr_m ) delete ptr_m;
}
void ManifoldRattleIntegrator::init()
{
// Makes sure the manifold params are set initially.
update_var_params();
reset_dt();
}
void ManifoldRattleIntegrator::update_params()
{
update_var_params();
}
void ManifoldRattleIntegrator::update_var_params()
{
if( nevery > 0 ){
stats.x_iters = 0;
stats.v_iters = 0;
stats.natoms = 0;
stats.x_iters_per_atom = 0.0;
stats.v_iters_per_atom = 0.0;
}
double **ptr_params = ptr_m->get_params();
for( int i = 0; i < nvars; ++i ){
if( is_var[i] ){
tvars[i] = input->variable->find(tstrs[i]);
if( tvars[i] < 0 ){
error->all(FLERR,
"Variable name for fix nve/manifold/rattle does not exist");
}
if( input->variable->equalstyle(tvars[i]) ){
tstyle[i] = EQUAL;
double new_val = input->variable->compute_equal(tvars[i]);
// fprintf( stdout, "New value of var %d is now %f\n", i+1, new_val );
*(ptr_params[i]) = new_val;
}else{
error->all(FLERR,
"Variable for fix nve/manifold/rattle is invalid style");
}
}
}
}
double ManifoldRattleIntegrator::memory_usage()
{
double bytes = 0.0;
bytes += sizeof(statistics);
bytes += sizeof(*ptr_m) + sizeof(ptr_m);
bytes += nvars*sizeof(double) + sizeof(double*);
bytes += nvars*( sizeof(char*) + 3*sizeof(int) );
return bytes;
}

View File

@ -0,0 +1,131 @@
/* ----------------------------------------------------------------------
Lammps - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
-----------------------------------------------------------------------
This file is a part of the USER-MANIFOLD package.
Copyright (2013-2015) Stefan Paquay, Eindhoven University of Technology.
License: GNU General Public License v.2
See the README file in the top-level LAMMPS directory.
This file is part of the user-manifold package written by
Stefan Paquay at the Eindhoven University of Technology.
This module makes it possible to do MD with particles constrained
to pretty arbitrary manifolds characterised by some constraint function
g(x,y,z) = 0 and its normal grad(g). The number of manifolds available
right now is limited but can be extended straightforwardly by making
a new class that inherits from manifold and implements all pure virtual
methods.
Thanks to Remy Kusters for beta-testing!
------------------------------------------------------------------------- */
#ifndef LMP_MANIFOLD_RATTLE_INTEGRATOR_H
#define LMP_MANIFOLD_RATTLE_INTEGRATOR_H
#include "manifold.h"
#include "pointers.h"
namespace LAMMPS_NS {
// A general util function for handling args.
void shift_args_back( int argi, char ***arg_ptr, int delta, int &narg );
class ManifoldRattleIntegratorBase : public Pointers
{
public:
struct statistics {
statistics() : x_iters(0), v_iters(0), x_iters_per_atom(0),
v_iters_per_atom(0), natoms(0), dofs_removed(0),
last_out(0) {}
double x_iters, v_iters;
double x_iters_per_atom;
double v_iters_per_atom;
int natoms;
int dofs_removed;
bigint last_out;
};
ManifoldRattleIntegratorBase( LAMMPS*, int &, char *** );
~ManifoldRattleIntegratorBase();
virtual void init() = 0;
virtual void update_params() = 0;
virtual double memory_usage() = 0;
virtual void nve_x_rattle(int, int);
virtual void nve_v_rattle(int, int);
virtual int nparams() = 0;
void init_manifold( int, char **, manifold * );
void print_stats( const char * );
void reset_dt();
int parse_arg( int , int &, char *** );
int is_legal_keyword( const char * );
int dof(int,int);
int nevery;
protected:
double dtv, dtf;
double tolerance;
int max_iter;
char **tstrs;
int *tvars;
int *tstyle;
int *is_var;
statistics stats;
int update_style;
int nvars;
private:
virtual void rattle_manifold_x( double*, double*, double*, double, double, int = 0 ) = 0;
virtual void rattle_manifold_v( double*, double*, double*, double, int = 0 ) = 0;
};
class ManifoldRattleIntegrator : public ManifoldRattleIntegratorBase
{
public:
ManifoldRattleIntegrator(LAMMPS *, int&, char *** );
virtual ~ManifoldRattleIntegrator();
virtual void init();
virtual void update_params();
virtual double memory_usage();
virtual int nparams(){ return ptr_m->nparams(); }
private:
void update_var_params();
virtual void rattle_manifold_x( double*, double*, double*, double, double, int = 0 );
virtual void rattle_manifold_v( double*, double*, double*, double, int = 0);
manifold *ptr_m;
};
}
#endif // LMP_MANIFOLD_RATTLE_INTEGRATOR_H

View File

@ -0,0 +1,54 @@
#ifndef LMP_MANIFOLD_SPHERE_H
#define LMP_MANIFOLD_SPHERE_H
#include "manifold.h"
namespace LAMMPS_NS {
// A sphere:
class manifold_sphere : public manifold, protected Pointers {
public:
enum { NPARAMS = 1 };
manifold_sphere( LAMMPS *lmp, int, char ** ) : Pointers(lmp){}
virtual ~manifold_sphere(){}
virtual double g( const double *x )
{
double R = params[0];
double r2 = x[0]*x[0] + x[1]*x[1] + x[2]*x[2];
return r2 - R*R;
}
virtual double g_and_n( const double *x, double *nn )
{
double R = params[0];
double r2 = x[0]*x[0] + x[1]*x[1] + x[2]*x[2];
nn[0] = 2*x[0];
nn[1] = 2*x[1];
nn[2] = 2*x[2];
return r2 - R*R;
}
virtual void n( const double *x, double *nn )
{
nn[0] = 2*x[0];
nn[1] = 2*x[1];
nn[2] = 2*x[2];
}
virtual void H( double *x, double h[3][3] )
{
h[0][1] = h[0][2] = h[1][0] = h[1][2] = h[2][0] = h[2][1] = 0.0;
h[0][0] = h[1][1] = h[2][2] = 2.0;
}
static const char* type(){ return "sphere"; }
virtual const char *id(){ return type(); }
static int expected_argc(){ return NPARAMS; }
virtual int nparams(){ return NPARAMS; }
};
}
#endif // LMP_MANIFOLD_SPHERE_H

View File

@ -0,0 +1,54 @@
#include "manifold_spine.h"
#include <math.h>
using namespace LAMMPS_NS;
manifold_spine::manifold_spine( LAMMPS *lmp, int argc, char **argv ) : Pointers(lmp)
{}
/*
* Equation for spine is:
*
* -x^2 - y^2 + (a^2 - (z/c)^2)*( 1 + A*sin(const(z) *z^2) )^4
* const(z) = (z < 0) ? B2 : B
* params[5] = { a, A, B, B2, c }
*/
double manifold_spine::g( const double *x )
{
double a = params[0];
double A = params[1];
double B = params[2];
double B2 = params[3];
double c = params[4];
if (x[2]>0){
return -1.0*(x[0]*x[0]+x[1]*x[1])+(a*a-x[2]*x[2]/(c*c))*(1.0+pow(A*sin(B *x[2]*x[2]),4));
}else{
return -1.0*(x[0]*x[0]+x[1]*x[1])+(a*a-x[2]*x[2]) *(1.0+pow(A*sin(B2*x[2]*x[2]),4));
}
}
void manifold_spine::n( const double *x, double *nn )
{
double a = params[0];
double A = params[1];
double B = params[2];
double B2 = params[3];
double c = params[4];
if (x[2]>0){
nn[0] = -2.0*x[0];
nn[1] = -2.0*x[1];
nn[2] = 8.0*A*A*A*A*B*x[2]*(a*a-(x[2]*x[2]/(c*c)))*cos(B*x[2]*x[2]) *
pow(sin(B*x[2]*x[2]),3)-2*x[2]*(1+pow(A*sin(B*x[2]*x[2]),4))/(c*c);
}else{
nn[0] = -2.0*x[0];
nn[1] = -2.0*x[1];
nn[2] = 8.0*A*A*A*A*B2*x[2]*(a*a-x[2]*x[2])*cos(B2*x[2]*x[2]) *
pow(sin(B2*x[2]*x[2]),3)-2*x[2]*(1+pow(A*sin(B2*x[2]*x[2]),4));
}
}

View File

@ -0,0 +1,28 @@
#ifndef LMP_MANIFOLD_SPINE_H
#define LMP_MANIFOLD_SPINE_H
#include "manifold.h"
namespace LAMMPS_NS {
// A dendritic spine approximation:
class manifold_spine : public manifold, protected Pointers {
public:
enum { NPARAMS = 5 }; // Number of parameters.
manifold_spine( LAMMPS *lmp, int, char ** );
virtual ~manifold_spine(){}
virtual double g ( const double *x );
virtual void n ( const double *x, double *nn );
static const char* type(){ return "spine"; }
virtual const char *id(){ return type(); }
static int expected_argc(){ return NPARAMS; }
virtual int nparams(){ return NPARAMS; }
};
}
#endif // LMP_MANIFOLD_SPINE_H

View File

@ -0,0 +1,171 @@
#include "comm.h"
#include "manifold_spine_opt.h"
#include "math_const.h"
#include <math.h>
using namespace LAMMPS_NS;
using namespace MathConst;
manifold_spine_opt::manifold_spine_opt( LAMMPS *lmp, int argc, char **argv ) : Pointers(lmp), lu_table(2048)
{
}
manifold_spine_opt::~manifold_spine_opt()
{
}
/*
* Equation for spine is:
*
* -x^2 - y^2 + (a^2 - (z/c)^2)*( 1 + A*sin(const(z) *z^2) )^4
* const(z) = (z < 0) ? B2 : B
* params[5] = { a, A, B, B2, c }
*
* These functions are however optimized quite a bit (but still slow =()
*/
double manifold_spine_opt::g( const double *x )
{
const double B = params[2];
const double B2 = params[3];
const double x2 = x[0]*x[0];
const double y2 = x[1]*x[1];
const double z2 = x[2]*x[2];
// Delay branching up to here, compute everything that is used by both before.
if( x[2] > 0 ){
const double z2_c2 = inv_c2 * z2;
const double Bz2 = B * z2;
const double s_Bz2 = lu_table.sin(Bz2);
const double s2_Bz2 = s_Bz2*s_Bz2;
const double s4_Bz2 = s2_Bz2*s2_Bz2;
const double one_p_A4s4 = 1.0 + A4*s4_Bz2;
const double a2_m_z2_c2 = a2 - z2_c2;
return -x2 - y2 + a2_m_z2_c2 * one_p_A4s4;
}else{
const double Bz2 = B2 * z2;
const double s_Bz2 = lu_table.sin(Bz2);
const double s2_Bz2 = s_Bz2*s_Bz2;
const double s3_Bz2 = s_Bz2*s2_Bz2;
const double s4_Bz2 = s2_Bz2*s2_Bz2;
const double one_p_A4s4 = (1.0 + A4*s4_Bz2);
const double a2_m_z2 = a2 - z2;
return -x2 - y2 + a2_m_z2 * one_p_A4s4;
}
}
void manifold_spine_opt::n( const double *x, double *nn )
{
const double B = params[2];
const double B2 = params[3];
nn[0] = -2*x[0];
nn[1] = -2*x[1];
const double z2 = x[2]*x[2];
// Delay branching up to here, compute everything that is used by both before.
if( x[2] > 0 ){
const double z2_c2 = inv_c2 * z2;
const double z_c2 = inv_c2 * x[2];
const double Bz2 = B * z2;
const double c_Bz2 = lu_table.cos(Bz2);
const double s_Bz2 = lu_table.sin(Bz2);
const double s2_Bz2 = s_Bz2*s_Bz2;
const double s3_Bz2 = s_Bz2*s2_Bz2;
const double s4_Bz2 = s2_Bz2*s2_Bz2;
const double one_p_A4s4 = 1.0 + A4*s4_Bz2;
const double a2_m_z2_c2 = a2 - z2_c2;
nn[2] = 8*A4*B*x[2] * a2_m_z2_c2 * c_Bz2 * s3_Bz2 - 2.0*z_c2 * one_p_A4s4;
}else{
const double Bz2 = B2 * z2;
const double c_Bz2 = lu_table.cos(Bz2);
const double s_Bz2 = lu_table.sin(Bz2);
const double s2_Bz2 = s_Bz2*s_Bz2;
const double s3_Bz2 = s_Bz2*s2_Bz2;
const double s4_Bz2 = s2_Bz2*s2_Bz2;
const double one_p_A4s4 = (1.0 + A4*s4_Bz2);
const double a2_m_z2 = a2 - z2;
nn[2] = 8*A4*B2*x[2] * a2_m_z2 * c_Bz2 * s3_Bz2 - 2.0*x[2] * one_p_A4s4;
}
}
double manifold_spine_opt::g_and_n( const double *x, double *nn )
{
const double B = params[2];
const double B2 = params[3];
nn[0] = -2*x[0];
nn[1] = -2*x[1];
const double x2 = x[0]*x[0];
const double y2 = x[1]*x[1];
const double z2 = x[2]*x[2];
// Delay branching up to here, compute everything that is used by both before.
if( x[2] > 0 ){
const double z2_c2 = inv_c2 * z2;
const double z_c2 = inv_c2 * x[2];
const double Bz2 = B * z2;
const double c_Bz2 = lu_table.cos(Bz2);
const double s_Bz2 = lu_table.sin(Bz2);
const double s2_Bz2 = s_Bz2*s_Bz2;
const double s3_Bz2 = s_Bz2*s2_Bz2;
const double s4_Bz2 = s2_Bz2*s2_Bz2;
const double one_p_A4s4 = 1.0 + A4*s4_Bz2;
const double a2_m_z2_c2 = a2 - z2_c2;
nn[2] = 8*A4*B*x[2] * a2_m_z2_c2 * c_Bz2 * s3_Bz2 - 2.0*z_c2 * one_p_A4s4;
return -x2 - y2 + a2_m_z2_c2 * one_p_A4s4;
}else{
const double Bz2 = B2 * z2;
const double c_Bz2 = lu_table.cos(Bz2);
const double s_Bz2 = lu_table.sin(Bz2);
const double s2_Bz2 = s_Bz2*s_Bz2;
const double s3_Bz2 = s_Bz2*s2_Bz2;
const double s4_Bz2 = s2_Bz2*s2_Bz2;
const double one_p_A4s4 = (1.0 + A4*s4_Bz2);
const double a2_m_z2 = a2 - z2;
nn[2] = 8*A4*B2*x[2] * a2_m_z2 * c_Bz2 * s3_Bz2 - 2.0*x[2] * one_p_A4s4;
return -x2 - y2 + a2_m_z2 * one_p_A4s4;
}
}
void manifold_spine_opt::post_param_init()
{
const double a = params[0];
const double A = params[1];
const double c = params[4];
c2 = c*c;
inv_c2 = 1.0 / c2;
a2 = a*a;
const double A2 = A*A;
A4 = A2*A2;
}

View File

@ -0,0 +1,35 @@
#ifndef LMP_MANIFOLD_SPINE_OPT_H
#define LMP_MANIFOLD_SPINE_OPT_H
#include "manifold.h"
#include "sincos_lookup.h"
namespace LAMMPS_NS {
// A dendritic spine approximation:
class manifold_spine_opt : public manifold, protected Pointers {
public:
enum { NPARAMS = 5 }; // Number of parameters.
manifold_spine_opt( LAMMPS *lmp, int, char ** );
virtual ~manifold_spine_opt();
virtual double g ( const double *x );
virtual void n ( const double *x, double *nn );
virtual double g_and_n( const double *x, double *nn );
static const char* type(){ return "spine/opt"; }
virtual const char *id(){ return type(); }
static int expected_argc(){ return NPARAMS; }
virtual int nparams(){ return NPARAMS; }
virtual void post_param_init();
private:
double c2, inv_c2, a2, A4;
sincos_lookup lu_table;
};
}
#endif // LMP_MANIFOLD_SPINE_OPT_H

View File

@ -0,0 +1,54 @@
#ifndef LMP_MANIFOLD_SUPERSPHERE_H
#define LMP_MANIFOLD_SUPERSPHERE_H
#include "manifold.h"
namespace LAMMPS_NS {
// A sphere:
class manifold_supersphere : public manifold, protected Pointers {
public:
enum { NPARAMS = 2 };
manifold_supersphere( LAMMPS *lmp, int, char ** ) : Pointers(lmp){}
virtual ~manifold_supersphere(){}
double my_sign( double a )
{
return (a > 0) - (a < 0);
}
virtual double g( const double *x )
{
double R = params[0];
double q = params[1];
double xx = fabs(x[0]);
double yy = fabs(x[1]);
double zz = fabs(x[2]);
double rr = pow(xx,q) + pow(yy,q) + pow(zz,q);
return rr - pow(R,q);
}
virtual void n( const double *x, double *nn )
{
double q = params[1];
double xx = fabs(x[0]);
double yy = fabs(x[1]);
double zz = fabs(x[2]);
nn[0] = q * my_sign(x[0])*pow(xx,q-1);
nn[1] = q * my_sign(x[1])*pow(yy,q-1);
nn[2] = q * my_sign(x[2])*pow(zz,q-1);
}
static const char* type(){ return "supersphere"; }
virtual const char *id(){ return type(); }
static int expected_argc(){ return NPARAMS; }
virtual int nparams(){ return NPARAMS; }
};
}
#endif

View File

@ -0,0 +1,382 @@
/* ----------------------------------------------------------------------
Lammps - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
-----------------------------------------------------------------------
This file is a part of the USER-MANIFOLD package.
Copyright (2013-2014) Stefan Paquay, Eindhoven University of Technology.
License: GNU General Public License v.2.
See the README file in the top-level LAMMPS directory.
This file is part of the user-manifold package written by
Stefan Paquay at the Eindhoven University of Technology.
This module makes it possible to do MD with particles constrained
to pretty arbitrary manifolds characterised by some constraint function
g(x,y,z) = 0 and its normal grad(g). The number of manifolds available
right now is limited but can be extended straightforwardly by making
a new class that inherits from manifold and implements all pure virtual
methods.
Thanks to Remy Kusters for beta-testing!
------------------------------------------------------------------------- */
#include "atom.h"
#include "comm.h"
#include "compute.h"
#include "error.h"
#include "fix.h"
#include "force.h"
#include "group.h"
#include "math.h"
#include "modify.h"
#include "neighbor.h"
#include "update.h"
#include "manifold_thermostat_nh.h"
#include "manifold_rattle_integrator.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum {CONSTANT,EQUAL};
enum {NOBIAS,BIAS};
ManifoldThermostatNH::ManifoldThermostatNH( LAMMPS *lmp, int &narg,
char ***arg_ptr, const char *f_id,
int argi, int igroup, int groupbit )
: Pointers(lmp), id_temp(NULL), temperature(NULL), eta(NULL), eta_dot(NULL),
eta_dotdot(NULL), eta_mass(NULL), fix_id(f_id),
igroup(igroup), groupbit(groupbit)
{
dthalf = dt4 = dt8 = 0;
t_start = t_stop = t_period = t_current = t_target = ke_target = 0.0;
t_freq = drag = tdrag_factor = 0;
boltz = force->boltz, nktv2p = force->nktv2p;
tdof = 0;
mtchain = 3;
factor_eta = 0.0;
which = got_temp = 0;
// argi is passed from somewhere else.
while( argi < narg ){
char **arg = *arg_ptr;
if( comm->me == 0 ){
fprintf(screen,"Parsing arg \"%s\"\n",arg[argi]);
}
if( is_legal_keyword( arg[argi] ) ){
int status = parse_arg(argi, narg, arg_ptr);
if( status ){
char msg[256];
sprintf(msg,"Error parsing arg \"%s\".\n", arg[argi]);
error->all(FLERR, msg);
}
}else{
++argi;
}
}
if( !got_temp ) error->all(FLERR,"Fix nvt/manifold/rattle needs 'temp'!");
if( t_period < 0.0 ){
error->all(FLERR,"Fix nvt/manifold/rattle damping parameter must be > 0.0");
}
// Create temperature compute:
int n = strlen(fix_id)+6;
id_temp = new char[n];
strcpy(id_temp,fix_id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = group->names[igroup];
newarg[2] = (char*) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
int icompute = modify->find_compute(id_temp);
if( icompute < 0 ){
error->all(FLERR,"Temperature ID for fix nvt/manifold/rattle "
"does not exist");
}
temperature = modify->compute[icompute];
if( temperature->tempbias ) which = BIAS;
else which = NOBIAS;
// Set t_freq from t_period
t_freq = 1.0 / t_period;
// Init Nosé-Hoover chain:
eta = new double[mtchain];
eta_dot = new double[mtchain+1];
eta_dotdot = new double[mtchain];
eta_mass = new double[mtchain];
eta_dot[mtchain] = 0.0;
eta_dot[mtchain] = 0.0;
for( int ich = 0; ich < mtchain; ++ich ){
eta[ich] = eta_dot[ich] = eta_dotdot[ich] = 0.0;
}
}
ManifoldThermostatNH::~ManifoldThermostatNH()
{
// Deallocate heap-allocated objects.
if( eta ) delete[] eta;
if( eta_dot ) delete[] eta_dot;
if( eta_dotdot ) delete[] eta_dotdot;
if( eta_mass ) delete[] eta_mass;
modify->delete_compute(id_temp);
if( id_temp ) delete[] id_temp;
}
int ManifoldThermostatNH::parse_arg( int argi, int &narg, char ***arg_ptr )
{
if( strcmp( (*arg_ptr)[argi], "temp") == 0 ){
if( argi+3 >= narg )
error->all(FLERR,"Keyword 'temp' needs 3 arguments");
t_start = force->numeric(FLERR,(*arg_ptr)[argi+1]);
t_stop = force->numeric(FLERR,(*arg_ptr)[argi+2]);
t_period = force->numeric(FLERR,(*arg_ptr)[argi+3]);
t_target = t_start;
got_temp = 1;
shift_args_back( argi, arg_ptr, 4, narg );
return 0;
}else if( strcmp( (*arg_ptr)[argi], "tchain" ) == 0 ){
if( argi+1 >= narg )
error->all(FLERR,"Keyword 'tchain' needs 1 argument");
mtchain = force->inumeric(FLERR, (*arg_ptr)[argi+1]);
shift_args_back( argi, arg_ptr, 2, narg );
return 0;
}
return -1;
}
int ManifoldThermostatNH::is_legal_keyword( const char *arg )
{
return (strcmp(arg, "temp") == 0) || (strcmp(arg, "tchain") == 0 );
}
void ManifoldThermostatNH::init()
{
int icompute = modify->find_compute(id_temp);
if( icompute < 0 ){
error->all(FLERR,"Temperature ID for fix nvt/manifold/rattle "
"does not exist");
}
temperature = modify->compute[icompute];
if( temperature->tempbias ) which = BIAS;
else which = NOBIAS;
}
double ManifoldThermostatNH::memory_usage()
{
double bytes = 0.0;
return bytes;
}
// All of this is copied/stolen from fix_nh:
void ManifoldThermostatNH::setup(int vflag)
{
compute_temp_target();
t_current = temperature->compute_scalar();
tdof = temperature->dof;
// Compute/set eta-masses:
double inv_t_freq2 = 1.0 / (t_freq*t_freq);
eta_mass[0] = tdof * boltz * t_target * inv_t_freq2;
for( int ich = 1; ich < mtchain; ++ich ){
eta_mass[ich] = boltz * t_target * inv_t_freq2;
}
for( int ich = 1; ich < mtchain; ++ich ){
eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1] -
boltz * t_target ) / eta_mass[ich];
}
}
void ManifoldThermostatNH::compute_temp_target()
{
t_current = temperature->compute_scalar();
tdof = temperature->dof;
double delta = update->ntimestep - update->beginstep;
if (delta != 0.0){
delta /= update->endstep - update->beginstep;
}
tdof = temperature->dof;
t_target = t_start + delta * (t_stop-t_start);
ke_target = tdof * boltz * t_target;
}
void ManifoldThermostatNH::nhc_temp_integrate()
{
int ich;
// t_current = temperature->compute_scalar();
// tdof = temperature->dof;
compute_temp_target();
double expfac, kecurrent = tdof * boltz * t_current;
double inv_t_freq2 = 1.0 / (t_freq*t_freq);
eta_mass[0] = tdof * boltz * t_target * inv_t_freq2;
for( int ich = 1; ich < mtchain; ++ich ){
eta_mass[ich] = boltz * t_target * inv_t_freq2;
}
if( eta_mass[0] > 0.0 ){
eta_dotdot[0] = (kecurrent - ke_target)/eta_mass[0];
}else{
eta_dotdot[0] = 0;
}
for( ich = mtchain-1; ich > 0; --ich ){
expfac = exp(-dt8*eta_dot[ich+1]);
eta_dot[ich] *= expfac;
eta_dot[ich] += eta_dotdot[ich] * dt4;
eta_dot[ich] *= tdrag_factor * expfac;
}
expfac = exp(-dt8*eta_dot[1]);
eta_dot[0] *= expfac;
eta_dot[0] += eta_dotdot[0] * dt4;
eta_dot[0] *= tdrag_factor * expfac;
factor_eta = exp(-dthalf*eta_dot[0]);
if( factor_eta == 0 ){
char msg[2048];
sprintf(msg, "WTF, factor_eta is 0! dthalf = %f, eta_dot[0] = %f",
dthalf, eta_dot[0]);
error->all(FLERR,msg);
}
nh_v_temp();
t_current *= factor_eta*factor_eta;
kecurrent = tdof * boltz * t_current;
if( eta_mass[0] > 0.0 ){
eta_dotdot[0] = (kecurrent - ke_target) / eta_mass[0];
}else{
eta_dotdot[0] = 0.0;
}
for( int ich = 1; ich < mtchain; ++ich ){
eta[ich] += dthalf*eta_dot[ich];
}
eta_dot[0] *= expfac;
eta_dot[0] += eta_dotdot[0]*dt4;
eta_dot[0] *= expfac;
for( int ich = 1; ich < mtchain; ++ich ){
expfac = exp(-dt8*eta_dot[ich+1]);
eta_dot[ich] *= expfac;
eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1]
- boltz*t_target) / eta_mass[ich];
eta_dot[ich] *= eta_dotdot[ich] * dt4;
eta_dot[ich] *= expfac;
}
}
void ManifoldThermostatNH::nh_v_temp()
{
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if( igroup == atom->firstgroup) nlocal = atom->nfirst;
if( which == NOBIAS ){
for( int i = 0; i < nlocal; ++i ){
if( mask[i] & groupbit ){
v[i][0] *= factor_eta;
v[i][1] *= factor_eta;
v[i][2] *= factor_eta;
}
}
}else if( which == BIAS ){
for( int i = 0; i < nlocal; ++i ){
if( mask[i] & groupbit ){
temperature->remove_bias(i,v[i]);
v[i][0] *= factor_eta;
v[i][1] *= factor_eta;
v[i][2] *= factor_eta;
temperature->restore_bias(i,v[i]);
}
}
}
}
void ManifoldThermostatNH::print_stats( const char *header )
{
}
void ManifoldThermostatNH::reset_dt()
{
dthalf = 0.5 * update->dt;
dt4 = 0.25 * update->dt;
dt8 = 0.125 * update->dt;
tdrag_factor = 1.0 - (update->dt * t_freq * drag);
}
void ManifoldThermostatNH::print_stuff()
{
if( comm->me != 0 ) return;
fprintf(screen, "\n\n********* Thermostat settings ************\n");
fprintf(screen, "Target temperature is %f\n", t_target);
fprintf(screen, "\n\n******************************************\n\n\n");
}
void ManifoldThermostatNH::update_t_current()
{
if (which == BIAS && neighbor->ago == 0)
t_current = temperature->compute_scalar();
}

View File

@ -0,0 +1,91 @@
/* ----------------------------------------------------------------------
Lammps - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
-----------------------------------------------------------------------
This file is a part of the USER-MANIFOLD package.
Copyright (2013-2015) Stefan Paquay, Eindhoven University of Technology.
License: GNU General Public License v.2
See the README file in the top-level LAMMPS directory.
This file is part of the user-manifold package written by
Stefan Paquay at the Eindhoven University of Technology.
This module makes it possible to do MD with particles constrained
to pretty arbitrary manifolds characterised by some constraint function
g(x,y,z) = 0 and its normal grad(g). The number of manifolds available
right now is limited but can be extended straightforwardly by making
a new class that inherits from manifold and implements all pure virtual
methods.
Thanks to Remy Kusters for beta-testing!
------------------------------------------------------------------------- */
#ifndef LMP_MANIFOLD_THERMOSTAT_NH_H
#define LMP_MANIFOLD_THERMOSTAT_NH_H
#include "pointers.h"
namespace LAMMPS_NS {
class ManifoldThermostatNH : protected Pointers {
public:
ManifoldThermostatNH(LAMMPS *, int&, char ***, const char *,
int, int, int );
~ManifoldThermostatNH();
int parse_arg( int , int &, char *** );
int is_legal_keyword( const char * );
double memory_usage();
void print_stats( const char * );
void init();
void setup(int);
void compute_temp_target();
void nhc_temp_integrate();
void nh_v_temp();
void reset_dt();
int nevery;
void print_stuff();
void update_t_current();
private:
// Fix NVT needs more params for Nose-Hoover chain
double dthalf, dt4, dt8;
char *id_temp;
class Compute* temperature;
double t_start,t_stop, t_period;
double t_current,t_target,ke_target;
double t_freq, drag, tdrag_factor;
double boltz,nktv2p,tdof;
double *eta,*eta_dot;
double *eta_dotdot;
double *eta_mass;
int mtchain;
double factor_eta;
int which, got_temp;
const char *fix_id;
int igroup, groupbit;
};
}
#endif // LMP_MANIFOLD_THERMOSTAT_NH_H

View File

@ -0,0 +1,619 @@
#include "manifold_thylakoid.h"
#include "math.h"
#include "comm.h"
#include "domain.h" // For some checks regarding the simulation box.
#define MANIFOLD_THYLAKOID_DEBUG
using namespace LAMMPS_NS;
manifold_thylakoid::manifold_thylakoid( LAMMPS *lmp, int narg, char ** arg)
: Pointers(lmp)
{
// You can NOT depend on proper construction of the domains in
// the constructor, because the params are set by the function that
// calls the factory-construction method. Instead, the constructing
// fix should call post_param_init();
}
manifold_thylakoid::~manifold_thylakoid()
{
for( std::size_t i = 0; i < parts.size(); ++i ){
delete parts[i];
}
}
void manifold_thylakoid::post_param_init()
{
// Set coefficients:
lT = 3; // Radius of cylinder edges of stacks
LT = 15; // Size of faces of cylinder stacks
pad = 3; // Padding to prevent improper look-ups
wB = 3.0;
LB = 10.0;
lB = 3.0;
wB = params[0];
LB = params[1];
lB = params[2];
if( comm->me == 0 ){
fprintf(screen,"My params are now: lT = %f, LT = %f, pad = %f, "
"wB = %f, LB = %f, lB = %f\n", lT, LT, pad, wB, LB, lB );
fprintf(screen,"Calling init_domains() from post_param_init().\n");
}
init_domains();
checkup();
}
void manifold_thylakoid::checkup()
{
if( comm->me == 0 ){
fprintf(screen,"This is checkup of thylakoid %p\n", this);
fprintf(screen,"I have %ld parts. They are:\n", parts.size());
for( int i = 0; i < parts.size(); ++i ){
fprintf(screen, "[%f, %f] x [%f, %f] x [%f, %f]\n",
parts[i]->xlo, parts[i]->xhi,
parts[i]->ylo, parts[i]->yhi,
parts[i]->zlo, parts[i]->zhi );
}
fprintf(screen,"My params are:\n");
for( int i = 0; i < NPARAMS; ++i ){
fprintf(screen,"%f\n", params[i]);
}
}
}
double manifold_thylakoid::g( const double *x )
{
int err = 0;
std::size_t idx;
thyla_part *p = get_thyla_part(x,&err,&idx);
if(err){
char msg[2048];
sprintf(msg,"Error getting thyla_part for x = (%f, %f, %f)",x[0],x[1],x[2]);
error->one(FLERR,msg);
}
double con_val = p->g(x);
if( isfinite(con_val) ){
return con_val;
}else{
char msg[2048];
sprintf(msg,"Error, thyla_part of type %s returned %f as constraint val!",
p->type, con_val);
error->one(FLERR,msg);
return 0;
}
}
void manifold_thylakoid::n( const double *x, double *n )
{
int err = 0;
std::size_t idx;
thyla_part *p = get_thyla_part(x,&err,&idx);
if(err){
char msg[2048];
sprintf(msg,"Error getting thyla_part for x = (%f, %f, %f)",x[0],x[1],x[2]);
error->one(FLERR,msg);
}
p->n(x,n);
if( isfinite(n[0]) && isfinite(n[1]) && isfinite(n[2]) ){
return;
}else{
char msg[2048];
sprintf(msg,"Error, thyla_part of type %s returned (%f,%f,%f) as gradient!",
p->type, n[0], n[1], n[2]);
error->one(FLERR,msg);
}
}
thyla_part *manifold_thylakoid::get_thyla_part( const double *x, int *err_flag, std::size_t *idx )
{
for( std::size_t i = 0; i < parts.size(); ++i ){
thyla_part *p = parts[i];
if( is_in_domain(p,x) ){
if( idx != NULL ) *idx = i;
return p;
}
}
char msg[2048];
sprintf(msg,"Could not find thyla_part for x = (%f,%f,%f)", x[0],x[1],x[2]);
error->one(FLERR,msg);
return NULL;
}
void manifold_thylakoid::init_domains()
{
if( wB + 2*lB > LT ){
char msg[2048];
sprintf(msg,"LT = %f not large enough to accomodate bridge with "
"wB = %f and lB = %f! %f > %f\n", LT, wB, lB, wB + 2*lB, LT);
error->one(FLERR,msg);
}
// Determine some constant coordinates:
x0 = -( 0.5*LB + lB + lT + LT + lT + pad);
y0 = -( 0.5*LT + lT + pad );
z0 = -15;
#ifdef MANIFOLD_THYLAKOID_DEBUG
if( comm->me == 0 ){
fprintf(screen,"x0, y0, z0 = %f, %f, %f\n",x0,y0,z0);
}
#endif // MANIFOLD_THYLAKOID_DEBUG
#ifndef USE_PHONY_LAMMPS
if( x0 < domain->boxlo[0] ){
char msg[2048];
sprintf(msg,"Thylakoid expects xlo of at most %f, but found %f",
x0, domain->boxlo[0]);
error->one(FLERR,msg);
}
if( y0 < domain->boxlo[1] ){
char msg[2048];
sprintf(msg,"Thylakoid expects ylo of at most %f, but found %f",
y0, domain->boxlo[1]);
error->one(FLERR,msg);
}
#endif
// Add some padding to prevent improper lookups.
z0 -= pad;
x1 = -x0;
y1 = -y0;
z1 = -z0;
Lx = x1 - x0;
Ly = y1 - y0;
Lz = z1 - z0;
#ifndef USE_PHONY_LAMMPS
char msg[2048];
if(x1 > domain->boxhi[0]){
sprintf(msg,"Expected xhi larger than current box has: %f > %f",
x1, domain->boxhi[0]);
error->one(FLERR,msg);
}
if(y1 > domain->boxhi[1]){
sprintf(msg,"Expected yhi larger than current box has: %f > %f",
y1, domain->boxhi[1]);
error->one(FLERR,msg);
}
// if(z1 > domain->boxhi[2]){
// sprintf(msg,"Expected zhi larger than current box has: %f > %f",
// z1, domain->boxhi[2]);
// error->one(FLERR,msg);
// }
#endif
// Create and add the manifold parts to the array.
thyla_part *p;
// Determine coordinates of domain boundaries and centres of "mass":
thyla_part_geom cllb, cllt, clrb, clrt; // Left thylakoid cylinder parts
thyla_part_geom pll, plb, plt, plr; // Left thylakoid plane parts
thyla_part_geom crlb, crlt, crrb, crrt; // Right thylakoid cylinder parts
thyla_part_geom prl, prb, prt, prr; // Right thylakoid plane parts
thyla_part_geom bl, br, bc; // Bridge left, right connectors and cylinder.
// The bridge is three parts.
// 1. A connector between the right face of the left grana and a cylinder
// 2. A connector between the left face of the right grana and a cylinder
// 3. The aforementioned cylinder.
// 1.:
// Args: pt, X0, R0, R, x0, y0, z0, sign
double rB = 0.5*wB;
double Reff = rB + lB;
double safety_fac = 1.0;
bl.pt[0] = -0.5*LB;
bl.pt[1] = 0;
bl.pt[2] = 0;
bl.lo[0] = bl.pt[0] - safety_fac*lB;
bl.lo[1] = -(1.0 + safety_fac) * Reff;
bl.lo[2] = -(1.0 + safety_fac) * Reff;
bl.hi[0] = bl.pt[0];
bl.hi[1] = (1.0 + safety_fac) * Reff;
bl.hi[2] = (1.0 + safety_fac) * Reff;
// double X0, double R0, double R, double s,
#ifdef MANIFOLD_THYLAKOID_DEBUG
if( comm->me == 0 ){
fprintf(screen,"x0, r0, R = %f, %f, %f\n", bl.pt[0], rB, lB);
}
#endif // MANIFOLD_THYLAKOID_DEBUG
p = make_cyl_to_plane_part(bl.pt[0], rB, lB, -1, bl.pt);
set_domain(p, bl.lo, bl.hi);
parts.push_back(p);
// 2.:
br.pt[0] = 0.5*LB;
br.pt[1] = 0;
br.pt[2] = 0;
br.lo[0] = br.pt[0];
br.lo[1] = -(1.0 + safety_fac) * Reff;
br.lo[2] = -(1.0 + safety_fac) * Reff;
br.hi[0] = br.pt[0] + safety_fac*lB;
br.hi[1] = (1.0 + safety_fac) * Reff;
br.hi[2] = (1.0 + safety_fac) * Reff;
// double X0, double R0, double R, double s,
#ifdef MANIFOLD_THYLAKOID_DEBUG
if( comm->me == 0 ){
fprintf(screen,"x0, r0, R = %f, %f, %f\n", br.pt[0], rB, lB);
}
#endif // MANIFOLD_THYLAKOID_DEBUG
p = make_cyl_to_plane_part(br.pt[0], rB, lB, 1, br.pt);
set_domain(p, br.lo, br.hi);
parts.push_back(p);
// 3.:
// Cylinder in between:
bc.pt[0] = 0;
bc.pt[1] = 0;
bc.pt[2] = 0;
bc.lo[0] = bl.pt[0];
bc.lo[1] = -Reff;
bc.lo[2] = -Reff;
bc.hi[0] = br.pt[0];
bc.hi[1] = Reff;
bc.hi[2] = Reff;
p = make_cyl_part( 0, 1, 1, bc.pt, rB );
set_domain( p, bc.lo, bc.hi );
#ifdef MANIFOLD_THYLAKOID_DEBUG
if( comm->me == 0 ){
fprintf(screen,"Cylinder lives on [ %f x %f ] x [ %f x %f ] x [ %f x %f]\n",
bc.lo[0], bc.hi[0], bc.lo[1], bc.hi[1], bc.lo[2], bc.hi[2]);
}
#endif // MANIFOLD_THYLAKOID_DEBUG
parts.push_back(p);
// The stack on the left:
cllb.lo[0] = x0;
cllb.lo[1] = y0;
cllb.lo[2] = z0;
cllb.pt[0] = x0 + pad + lT;
cllb.pt[1] = y0 + pad + lT;
cllb.pt[2] = 0;
cllb.hi[0] = cllb.pt[0];
cllb.hi[1] = cllb.pt[1];
cllb.hi[2] = z1;
p = make_cyl_part(1,1,0,cllb.pt,lT);
set_domain(p, cllb.lo, cllb.hi);
parts.push_back(p);
// left left top cylinder
cllt = cllb;
cllt.lo[1] = y1 - pad - lT;
cllt.hi[1] = y1;
cllt.pt[1] = cllb.pt[1] + LT;
p = make_cyl_part(1,1,0,cllt.pt,lT);
set_domain(p, cllt.lo, cllt.hi);
parts.push_back(p);
// left right bottom cylinder
clrb = cllb;
clrb.pt[0] += LT;
clrb.lo[0] = clrb.pt[0];
clrb.hi[0] = clrb.lo[0] + lT + lB;
p = make_cyl_part(1,1,0,clrb.pt,lT);
set_domain(p, clrb.lo, clrb.hi);
parts.push_back(p);
// left right top cylinder
clrt = clrb;
clrt.pt[1] += LT;
clrt.lo[1] = y1 - pad - lT;
clrt.hi[1] = y1;
p = make_cyl_part(1,1,0,clrt.pt,lT);
set_domain(p, clrt.lo, clrt.hi);
parts.push_back(p);
// left left plane
pll.pt[0] = x0 + pad;
pll.pt[1] = 0;
pll.pt[2] = 0;
pll.lo[0] = x0;
pll.lo[1] = cllb.pt[1];
pll.lo[2] = z0;
pll.hi[0] = pll.lo[0] + pad + lT;
pll.hi[1] = pll.lo[1] + LT;
pll.hi[2] = z1;
p = make_plane_part(1,0,0,pll.pt);
set_domain(p, pll.lo, pll.hi);
parts.push_back(p);
// left bottom plane
plb.pt[0] = x0 + pad + lT + 0.5*LT;
plb.pt[1] = y0 + pad;
plb.pt[2] = 0;
plb.lo[0] = x0 + pad + lT;
plb.lo[1] = y0;
plb.lo[2] = z0;
plb.hi[0] = plb.lo[0] + LT;
plb.hi[1] = plb.lo[1] + pad + lT;
plb.hi[2] = z1;
p = make_plane_part(0,1,0,plb.pt);
set_domain(p, plb.lo, plb.hi);
parts.push_back(p);
// left top plane
plt = plb;
plt.lo[1] = cllb.pt[1] + LT;
plt.hi[1] = y1;
plt.pt[1] = y1 - pad;
p = make_plane_part(0,1,0,plt.pt);
set_domain(p, plt.lo, plt.hi);
parts.push_back(p);
// left right plane
plr = pll;
plr.lo[0] = bl.lo[0] - 0.5;
plr.lo[1] = y0 - pad;
plr.hi[0] = bl.lo[0] + 0.5;
plr.hi[1] = y1 + pad;
plr.pt[0] = bl.pt[0] - lB;
plr.pt[1] = 0.0;
plr.pt[2] = 0.0;
plr.hi[2] = z1 + pad;
plr.lo[2] = z0 - pad;
p = make_plane_part(1,0,0,plr.pt);
set_domain(p, plr.lo, plr.hi);
parts.push_back(p);
// Check if this plane lines up with bl:
if( fabs(plr.pt[0] - bl.pt[0] + lB) > 1e-8 ){
char msg[2048];
sprintf(msg,"Origins of plane left right and bridge left misaligned! %f != %f!\n",
plr.pt[0], bl.pt[0] - lB );
error->one(FLERR,msg);
}
// Now, for the right stack, you can mirror the other...
// To mirror them you need to invert lo[0] and hi[0] and flip their sign.
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &crlb, &cllb);
p = make_cyl_part(1,1,0,crlb.pt,lT);
set_domain(p, crlb.lo, crlb.hi);
parts.push_back(p);
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &crlt, &cllt);
p = make_cyl_part(1,1,0,crlt.pt,lT);
set_domain(p, crlt.lo, crlt.hi);
parts.push_back(p);
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &crrb, &clrb);
p = make_cyl_part(1,1,0,crrb.pt,lT);
set_domain(p, crrb.lo, crrb.hi);
parts.push_back(p);
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &crrt, &clrt);
p = make_cyl_part(1,1,0,crrt.pt,lT);
set_domain(p, crrt.lo, crrt.hi);
parts.push_back(p);
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &prl , &pll );
p = make_plane_part(1,0,0,prl.pt);
set_domain(p, prl.lo, prl.hi);
parts.push_back(p);
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &prb , &plb );
p = make_plane_part(0,1,0,prb.pt);
set_domain(p, prb.lo, prb.hi);
parts.push_back(p);
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &prt , &plt );
p = make_plane_part(0,1,0,prt.pt);
set_domain(p, prt.lo, prt.hi);
parts.push_back(p);
// Careful, this one is wrongly named.
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &prr, &plr);
p = make_plane_part(1,0,0,prr.pt);
set_domain(p, prr.lo, prr.hi);
parts.push_back(p);
if( fabs(prr.pt[0] - br.pt[0] - lB) > 1e-8 ){
char msg[2048];
sprintf(msg,"Origins of plane left right and bridge left misaligned! %f != %f!\n",
prr.pt[0], br.pt[0] + lB);
error->one(FLERR,msg);
}
// For debugging, print the domains and coms:
#ifdef MANIFOLD_THYLAKOID_DEBUG
if( comm->me == 0 ){
FILE *fp_doms = fopen("test_doms.dat","w");
FILE *fp_coms = fopen("test_coms.dat","w");
print_part_data(fp_doms, fp_coms);
fclose(fp_doms);
fclose(fp_coms);
}
#endif // MANIFOLD_THYLAKOID_DEBUG
}
void manifold_thylakoid::set_domain( thyla_part *p, const std::vector<double> &lo,
const std::vector<double> &hi )
{
#ifdef MANIFOLD_THYLAKOID_DEBUG
if( comm->me == 0 ){
fprintf(screen,"Adding part with domain [%f, %f] x [%f, %f] x [%f, %f]\n",
lo[0],hi[0],lo[1],hi[1],lo[2],hi[2] );
}
#endif // MANIFOLD_THYLAKOID_DEBUG
if( lo[0] >= hi[0] ){
char msg[2048];
sprintf(msg,"xlo >= xhi (%f >= %f)",lo[0],hi[0]);
error->one(FLERR,msg);
}else if( lo[1] >= hi[1] ){
char msg[2048];
sprintf(msg,"ylo >= yhi (%f >= %f)",lo[1],hi[1]);
error->one(FLERR,msg);
}else if( lo[2] >= hi[2] ){
char msg[2048];
sprintf(msg,"zlo >= zhi (%f >= %f)",lo[2],hi[2]);
error->one(FLERR,msg);
}
p->xlo = lo[0];
p->ylo = lo[1];
p->zlo = lo[2];
p->xhi = hi[0];
p->yhi = hi[1];
p->zhi = hi[2];
}
int manifold_thylakoid::is_in_domain( thyla_part *part, const double *x )
{
bool domain_ok = (x[0] >= part->xlo) && (x[0] <= part->xhi) &&
(x[1] >= part->ylo) && (x[1] <= part->yhi) &&
(x[2] >= part->zlo) && (x[2] <= part->zhi);
if( !domain_ok ) return false;
// From here on out, domain is ok.
if( part->type == thyla_part::THYLA_TYPE_CYL_TO_PLANE ){
double R0 = part->params[1];
double R = part->params[2];
double X0 = part->params[0];
double y = x[1];
double z = x[2];
double dist2 = y*y + z*z;
double RR = R0+R;
double RR2 = RR*RR;
double s = part->params[6];
if( dist2 < RR2 ){
return true;
}else{
// Domain was ok, but radius not.
return false;
}
}else{
return true;
}
}
thyla_part *manifold_thylakoid::make_plane_part (double a, double b, double c,
const std::vector<double> &pt )
{
double args[7];
args[0] = a;
args[1] = b;
args[2] = c;
args[3] = pt[0];
args[4] = pt[1];
args[5] = pt[2];
thyla_part *p = new thyla_part(thyla_part::THYLA_TYPE_PLANE,args,0,0,0,0,0,0);
return p;
}
thyla_part *manifold_thylakoid::make_cyl_part (double a, double b, double c,
const std::vector<double> &pt, double R)
{
double args[7];
args[0] = a;
args[1] = b;
args[2] = c;
args[3] = pt[0];
args[4] = pt[1];
args[5] = pt[2];
args[6] = R;
thyla_part *p = new thyla_part(thyla_part::THYLA_TYPE_CYL,args,0,0,0,0,0,0);
return p;
}
thyla_part *manifold_thylakoid::make_sphere_part(const std::vector<double> &pt, double R)
{
double args[7];
args[0] = R;
args[1] = pt[0];
args[2] = pt[1];
args[3] = pt[2];
thyla_part *p = new thyla_part(thyla_part::THYLA_TYPE_SPHERE,args,0,0,0,0,0,0);
return p;
}
thyla_part *manifold_thylakoid::make_cyl_to_plane_part(double X0, double R0, double R,
double s, const std::vector<double> &pt )
{
double args[7];
args[0] = X0;
args[1] = R0;
args[2] = R;
args[3] = pt[0];
args[4] = pt[1];
args[5] = pt[2];
args[6] = s;
thyla_part *p = new thyla_part(thyla_part::THYLA_TYPE_CYL_TO_PLANE,args,0,0,0,0,0,0);
return p;
}
void manifold_thylakoid::print_part_data( FILE *fp_doms, FILE *fp_coms )
{
for( std::size_t i = 0; i < parts.size(); ++i ){
thyla_part *p = parts[i];
fprintf(fp_doms, "%f %f\n", p->xlo, p->ylo);
fprintf(fp_doms, "%f %f\n", p->xlo, p->yhi);
fprintf(fp_doms, "%f %f\n", p->xhi, p->yhi);
fprintf(fp_doms, "%f %f\n", p->xhi, p->ylo);
fprintf(fp_doms, "%f %f\n\n",p->xlo, p->ylo);
fprintf(fp_coms, "%f %f\n", p->x0, p->y0 );
}
}

View File

@ -0,0 +1,101 @@
#ifndef LMP_MANIFOLD_THYLAKOID_H
#define LMP_MANIFOLD_THYLAKOID_H
#include "manifold.h"
#include <cassert>
#include <vector>
#include "manifold_thylakoid_shared.h"
namespace LAMMPS_NS {
using LAMMPS_NS::thyla_part;
using LAMMPS_NS::thyla_part_geom;
class manifold_thylakoid : public manifold, protected Pointers {
public:
enum { NPARAMS = 3 };
manifold_thylakoid( LAMMPS *lmp, int, char ** );
virtual ~manifold_thylakoid();
virtual double g( const double *x );
virtual void n( const double *x, double *n );
static const char* type(){ return "thylakoid"; }
virtual const char *id(){ return type(); }
static int expected_argc(){ return NPARAMS; }
virtual int nparams(){ return NPARAMS; }
virtual void post_param_init();
virtual void checkup(); // Some diagnostics...
private:
void init_domains();
thyla_part *get_thyla_part( const double *x, int *err_flag, std::size_t *idx = NULL );
int is_in_domain( thyla_part *p, const double *x );
void check_overlap();
std::vector<thyla_part*> parts;
thyla_part *make_plane_part (double a, double b, double c,
const std::vector<double> &pt);
thyla_part *make_cyl_part (double a, double b, double c,
const std::vector<double> &pt, double R);
thyla_part *make_sphere_part(const std::vector<double> &pt, double R);
thyla_part *make_cyl_to_plane_part(double X0, double R0, double R, double s,
const std::vector<double> &pt );
void set_domain( thyla_part *p, const std::vector<double> &lo,
const std::vector<double> &hi );
void add_part( thyla_part *part, double x0, double y0,
double z0, double x1, double y1, double z1 )
{
int me = -1;
MPI_Comm_rank(world, &me);
std::vector<double> lo(3), hi(3);
lo[0] = x0; lo[1] = y0; lo[2] = z0;
hi[0] = x1; hi[1] = y1; hi[2] = z1;
set_domain(part, lo, hi);
parts.push_back( part );
fprintf(my_screen,"This is proc %d: Added part %p of type %d, "
"size is now %lu.\n", me, part, part->type, parts.size());
fflush(my_screen);
}
void print_part_data( FILE *fp_doms, FILE *fp_coms );
void assert_part_ok( const thyla_part *part ){
if( !part ) error->one(FLERR,"Error! Part was NULL!");
}
// Coefficients for the thylakoid model. At the moment it is just
// a cylinder, we slowly expand it.
double pad; // Padding to make sure periodic images are mapped back properly.
FILE *my_screen;
double LB, lT, lB, wB, LT;
// Domain size:
double x0, y0, z0;
double x1, y1, z1;
double Lx, Ly, Lz;
};
} // namespace LAMMPS_NS
#endif // LMP_MANIFOLD_THYLAKOID_H

View File

@ -0,0 +1,234 @@
#include "manifold_thylakoid_shared.h"
using namespace LAMMPS_NS;
thyla_part::thyla_part( int type, double *args, double xlo, double ylo, double zlo,
double xhi, double yhi, double zhi )
: type(type), xlo(xlo), xhi(xhi),
ylo(ylo), yhi(yhi), zlo(zlo), zhi(zhi),
x0( (type == THYLA_TYPE_SPHERE) ? params[1] : params[3] ),
y0( (type == THYLA_TYPE_SPHERE) ? params[2] : params[4] ),
z0( (type == THYLA_TYPE_SPHERE) ? params[3] : params[5] )
{
switch(type){
case THYLA_TYPE_PLANE: // a*(x-x0) + b*(y-y0) + c*(z-z0) = 0
params[0] = args[0]; // a
params[1] = args[1]; // b
params[2] = args[2]; // c
params[3] = args[3]; // x0
params[4] = args[4]; // y0
params[5] = args[5]; // z0
break;
case THYLA_TYPE_SPHERE: // (x-x0)^2 + (y-y0)^2 + (z-z0)^2 - R^2 = 0
params[0] = args[0]; // R
params[1] = args[1]; // x0
params[2] = args[2]; // y0
params[3] = args[3]; // z0
break;
case THYLA_TYPE_CYL: // a*(x-x0)^2 + b*(y-y0)^2 + c*(z-z0)^2 - R^2 = 0
params[0] = args[0]; // a
params[1] = args[1]; // b
params[2] = args[2]; // c
params[3] = args[3]; // x0
params[4] = args[4]; // y0
params[5] = args[5]; // z0
params[6] = args[6]; // R
if( (args[0] != 0.0) && (args[1] != 0.0) && (args[2] != 0.0) ){
err_flag = -1;
return;
}
// The others should be 1.
if( (args[0] != 1.0) && (args[0] != 0.0) &&
(args[1] != 1.0) && (args[1] != 0.0) &&
(args[2] != 1.0) && (args[2] != 0.0) ){
err_flag = -1;
}
break;
case THYLA_TYPE_CYL_TO_PLANE:
/*
* Funky bit that connects a cylinder to a plane.
* It is what you get by rotating the equation
* r(x) = R0 + R*( 1 - sqrt( 1 - ( ( X0 - x ) /R )^2 ) ) around the x-axis.
* I kid you not.
*
* The shape is symmetric so you have to set whether it is the "left" or
* "right" end by truncating it properly. It is designed to run from
* X0 to X0 + R if "right" or X0 - R to X0 if "left".
*
* As params it takes X0, R0, R, and a sign that determines whether it is
* "left" (args[3] < 0) or "right" (args[3] > 0).
*
* The args are: X0, R0, R, x0, y0, z0, sign.
*/
params[0] = args[0];
params[1] = args[1];
params[2] = args[2];
params[3] = args[3];
params[4] = args[4];
params[5] = args[5];
params[6] = args[6];
break;
default:
err_flag = -1;
}
}
thyla_part::~thyla_part()
{}
double thyla_part::g(const double *x)
{
switch(type){
case THYLA_TYPE_PLANE:{ // a*(x-x0) + b*(y-y0) + c*(z-z0) = 0
double a = params[0];
double b = params[1];
double c = params[2];
double dx = x[0] - params[3];
double dy = x[1] - params[4];
double dz = x[2] - params[5];
return a*dx + b*dy + c*dz;
break;
}
case THYLA_TYPE_SPHERE:{ // (x-x0)^2 + (y-y0)^2 + (z-z0)^2 - R^2 = 0
double R2 = params[0]*params[0];
double dx = x[0] - params[1];
double dy = x[1] - params[2];
double dz = x[2] - params[3];
return dx*dx + dy*dy + dz*dz - R2;
break;
}
case THYLA_TYPE_CYL:{ // a*(x-x0)^2 + b*(y-y0)^2 + c*(z-z0)^2 - R^2 = 0
double a = params[0];
double b = params[1];
double c = params[2];
double X0 = params[3];
double Y0 = params[4];
double Z0 = params[5];
double R = params[6];
double dx = x[0] - X0;
double dy = x[1] - Y0;
double dz = x[2] - Z0;
return a*dx*dx + b*dy*dy + c*dz*dz - R*R;
break;
}
case THYLA_TYPE_CYL_TO_PLANE:{
double X0 = params[0];
double R0 = params[1];
double R = params[2];
double s = sign(params[6]);
// Determine the centre of the sphere.
double dx = (x[0] - X0);
double dyz = sqrt( x[1]*x[1] + x[2]*x[2] );
double rdyz = dyz - (R0 + R);
double norm = 1.0 / (2.0*R);
// Maybe sign is important here...
double g = norm*(dx*dx + rdyz*rdyz - R*R);
return g;
}
default:
err_flag = -1;
return 0; // Mostly to get rid of compiler werrors.
break;
}
}
void thyla_part::n( const double *x, double *n )
{
switch(type){
case THYLA_TYPE_PLANE:{ // a*(x-x0) + b*(y-y0) + c*(z-z0) = 0
double a = params[0];
double b = params[1];
double c = params[2];
n[0] = a;
n[1] = b;
n[2] = c;
break;
}
case THYLA_TYPE_SPHERE:{ // (x-x0)^2 + (y-y0)^2 + (z-z0)^2 - R^2 = 0
double dx = x[0] - params[1];
double dy = x[1] - params[2];
double dz = x[2] - params[3];
n[0] = 2*dx;
n[1] = 2*dy;
n[2] = 2*dz;
break;
}
case THYLA_TYPE_CYL:{ // a*(x-x0)^2 + b*(y-y0)^2 + c*(z-z0)^2 - R^2 = 0
double a = params[0];
double b = params[1];
double c = params[2];
double X0 = params[3];
double Y0 = params[4];
double Z0 = params[5];
double R = params[6];
double dx = x[0] - X0;
double dy = x[1] - Y0;
double dz = x[2] - Z0;
n[0] = 2*a*dx;
n[1] = 2*b*dy;
n[2] = 2*c*dz;
break;
}
case THYLA_TYPE_CYL_TO_PLANE:{
double X0 = params[0];
double R0 = params[1];
double R = params[2];
double s = sign(params[6]);
// Determine the centre of the sphere.
double dx = s*(x[0] - X0);
double ryz = sqrt( x[1]*x[1] + x[2]*x[2] );
// Maybe sign is important here...
// Normalize g and n so that the normal is continuous:
double norm = 1.0 / (2.0 * R);
double x0[3];
x0[0] = X0-R;
x0[1] = 0;
x0[2] = R0+R;
n[0] = s*2*dx*norm;
double const_part = 1.0 - (R0 + R) / ryz;
n[1] = 2*x[1]*const_part*norm;
n[2] = 2*x[2]*const_part*norm;
break;
}
default:
err_flag = -1;
break;
}
}
void thyla_part_geom::mirror( unsigned int axis, thyla_part_geom *m,
const thyla_part_geom *o )
{
// Since dir is already the index of the array this is really simple:
m->lo[0] = o->lo[0];
m->lo[1] = o->lo[1];
m->lo[2] = o->lo[2];
m->pt[0] = o->pt[0];
m->pt[1] = o->pt[1];
m->pt[2] = o->pt[2];
m->hi[0] = o->hi[0];
m->hi[1] = o->hi[1];
m->hi[2] = o->hi[2];
m->lo[axis] = -o->hi[axis];
m->hi[axis] = -o->lo[axis];
m->pt[axis] = -o->pt[axis];
}

View File

@ -0,0 +1,74 @@
#ifndef MANIFOLD_THYLAKOID_SHARED_H
#define MANIFOLD_THYLAKOID_SHARED_H
#include <vector>
#include <cassert>
#include "math.h"
#include "pointers.h"
#include "error.h"
#include "lammps.h"
#include "manifold.h"
namespace LAMMPS_NS {
// The thylakoid is composed of many parts
struct thyla_part {
enum thyla_types {
THYLA_TYPE_PLANE,
THYLA_TYPE_SPHERE,
THYLA_TYPE_CYL,
THYLA_TYPE_CYL_TO_PLANE
};
thyla_part( int type, double *args, double xlo, double ylo, double zlo,
double xhi, double yhi, double zhi );
thyla_part() : type(-1), x0(-1337), y0(-1337), z0(-1337){}
~thyla_part();
double g( const double *x );
void n( const double *x, double *n );
int type;
double params[7];
double tol;
int maxit;
int err_flag;
tagint a_id;
double xlo, xhi, ylo, yhi, zlo, zhi;
double x0, y0, z0;
}; // struct thyla_part
struct thyla_part_geom {
thyla_part_geom() : pt(3), lo(3), hi(3){}
std::vector<double> pt, lo, hi;
// Function for mirroring thyla_geoms:
enum DIRS { DIR_X, DIR_Y, DIR_Z };
static void mirror( unsigned int axis, thyla_part_geom *m,
const thyla_part_geom *o );
}; // struct thyla_part_geom
} // namespace LAMMPS_NS
inline double sign( const double x )
{
return (x > 0.0) ? 1.0 : -1.0;
}
#endif // MANIFOLD_THYLAKOID_SHARED_H

View File

@ -0,0 +1,37 @@
#include "manifold_torus.h"
#include <math.h>
using namespace LAMMPS_NS;
manifold_torus::manifold_torus( LAMMPS *lmp, int argc, char **argv ) : Pointers(lmp)
{}
double manifold_torus::g( const double *x )
{
double R = params[0];
double r = params[1];
if( R < r ){
error->all(FLERR,"Large radius < small radius!");
}
double rad = sqrt(x[0]*x[0] + x[1]*x[1]);
double c = R - rad;
return c*c + x[2]*x[2] - r*r;
}
void manifold_torus::n( const double *x, double *n )
{
double R = params[0];
double r = params[1];
if( R < r ){
error->all(FLERR,"Large radius < small radius!");
}
double rad = sqrt(x[0]*x[0] + x[1]*x[1]);
double c = R - rad;
double fac = c / rad;
n[0] = -2.0 * fac * x[0];
n[1] = -2.0 * fac * x[1];
n[2] = 2.0 * x[2];
}

View File

@ -0,0 +1,25 @@
#ifndef LMP_MANIFOLD_TORUS_H
#define LMP_MANIFOLD_TORUS_H
#include "manifold.h"
namespace LAMMPS_NS {
class manifold_torus : public manifold, protected Pointers {
public:
enum {NPARAMS=2};
manifold_torus( LAMMPS *, int, char ** );
~manifold_torus(){}
virtual double g( const double *x );
virtual void n( const double *x, double *n );
static const char *type(){ return "torus"; }
virtual const char *id(){ return type(); }
static int expected_argc(){ return NPARAMS; }
virtual int nparams(){ return NPARAMS; }
};
}
#endif // LMP_MANIFOLD_TORUS_H

View File

@ -0,0 +1,102 @@
#include "sincos_lookup.h"
#include <cmath>
#include <iostream>
#define SINCOS_LOOKUP_NO_SAFETY
#if __cplusplus < 201103L
#define nullptr NULL
#endif // __cplusplus
sincos_lookup::sincos_lookup( int N ) : cos_table(nullptr), size(N)
{
std::cout << "Constructing table of size " << size
<< ", this might take a while... ";
make_cos_table();
std::cout << "Done!\n";
}
sincos_lookup::~sincos_lookup()
{
if( cos_table ) delete [] cos_table;
}
double sincos_lookup::sin( double x )
{
return this->cos( my_pi2 - x );
}
double sincos_lookup::cos( double x )
{
double x0 = x;
// Wrap to first quadrant:
x = fabs(x);
const double rounded_n_over_2pi = static_cast<int>(x / my_2pi) * my_2pi;
x -= rounded_n_over_2pi;
if( x > my_pi ){
// make use of cos( x ) = -cos(pi-x)
x -= 2.0*( x - my_pi );
}
// Get index:
int idx = static_cast<int>( x / dx_tab );
#ifdef SINCOS_LOOKUP_NO_SAFETY
// Assume safe:
#else
if( idx > size || idx < 0 ){
std::cerr << "Illegal table lookup for index " << idx << ", x="
<< x << ", x0=" << x0 << "!\n";
}
#endif // SINCOS_LOOKUP_NO_SAFETY
// static_cast always truncates to left, so lookup with "delta to right":
double delta = x - idx*dx_tab;
double cos_idx = cos_table[idx];
double diff_cos = cos_table[idx+1] - cos_idx;
return cos_idx + delta*diff_cos*inv_dx_tab;
}
void sincos_lookup::sincos(double x, double *s, double *c)
{
*s = this->sin(x);
*c = this->cos(x);
}
void sincos_lookup::test_lookup( )
{
double dx = 1e-5;
double max_err_c = 0, max_err_s = 0;
for( double x = 0; x < 2*my_2pi; x += dx ){
double c = this->cos(x);
double s = this->sin(x);
// Compare with math.h:
double c_err = c - cos(x);
double s_err = s - sin(x);
if( fabs(max_err_c) < fabs(c_err) ) max_err_c = c_err;
if( fabs(max_err_s) < fabs(s_err) ) max_err_s = s_err;
}
std::cout << "Largest lookup table errors (sin/cos): ( " << max_err_s
<< ", " << max_err_c << " ).\n";
}
void sincos_lookup::make_cos_table()
{
cos_table = new double [size+1];
if( !cos_table ){
// error.
}
dx_tab = my_pi/static_cast<double>(size-1);
inv_dx_tab = 1.0 / dx_tab;
for( int i = 0; i <= size; ++i ){
double f = i*dx_tab;
cos_table[i] = std::cos(f);
}
}

View File

@ -0,0 +1,35 @@
#ifndef SINCOS_LOOKUP_H
#define SINCOS_LOOKUP_H
class sincos_lookup
{
public:
sincos_lookup( int N );
~sincos_lookup();
double sin(double x);
double cos(double x);
void sincos(double x, double *s, double *c);
void test_lookup();
#if __cplusplus < 201103L
#define my_pi2 1.57079632679489661923
#define my_pi 3.14159265358979323846
#define my_2pi 6.28318530717958647692
#else
static constexpr const double my_pi2 = 1.57079632679489661923;
static constexpr const double my_pi = 3.14159265358979323846;
static constexpr const double my_2pi = 6.28318530717958647692;
#endif
private:
double *cos_table;
double dx_tab, inv_dx_tab;
int size;
void make_cos_table();
};
#endif // SINCOS_LOOKUP_H