Merge pull request #283 from akohlmey/grem-feature

gREM generalized replica exchange feature for USER-MISC
This commit is contained in:
sjplimp 2016-11-22 08:15:35 -07:00 committed by GitHub
commit d7b542101a
43 changed files with 17902 additions and 11 deletions

BIN
doc/src/Eqs/fix_grem.jpg Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 6.1 KiB

9
doc/src/Eqs/fix_grem.tex Normal file
View File

@ -0,0 +1,9 @@
\documentclass[12pt]{article}
\begin{document}
$$
T_{eff} = \lambda + \eta (H - H_0)
$$
\end{document}

View File

@ -687,6 +687,7 @@ package"_Section_start.html#start_3.
"eos/table/rx"_fix_eos_table_rx.html,
"flow/gauss"_fix_flow_gauss.html,
"gle"_fix_gle.html,
"grem"_fix_grem.html,
"imd"_fix_imd.html,
"ipi"_fix_ipi.html,
"langevin/drude"_fix_langevin_drude.html,

111
doc/src/fix_grem.txt Normal file
View File

@ -0,0 +1,111 @@
"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 grem command :h3
[Syntax:]
fix ID group-ID grem lambda eta H0 thermostat-ID :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
grem = style name of this fix command :l
lambda = intercept parameter of linear effective temperature function :l
eta = slope parameter of linear effective temperature function :l
H0 = shift parameter of linear effective temperature function :l
thermostat-ID = ID of Nose-Hoover thermostat or barostat used in simulation :l,ule
[Examples:]
fix fxgREM all grem 400 -0.01 -30000 fxnpt
thermo_modify press fxgREM_press :pre
fix fxgREM all grem 502 -0.15 -80000 fxnvt :pre
[Description:]
This fix implements the molecular dynamics version of the generalized
replica exchange method (gREM) originally developed by "(Kim)"_#Kim,
which uses non-Boltzmann ensembles to sample over first order phase
transitions. The is done by defining replicas with an enthalpy dependent
effective temperature
:c,image(Eqs/fix_grem.jpg)
with {eta} negative and steep enough to only intersect the characteristic
microcanonical temperature (Ts) of the system once, ensuring a unimodal
enthalpy distribution in that replica. {Lambda} is the intercept and
effects the generalized ensemble similar to how temperature effects a
Boltzmann ensemble. {H0} is a reference enthalpy, and is typically set
as the lowest desired sampled enthalpy. Further explanation can be
found in our recent papers "(Malolepsza)"_#Malolepsza.
This fix requires a Nose-Hoover thermostat fix reference passed to the
grem as {thermostat-ID}. Two distinct temperatures exist in this generalized
ensemble, the effective temperature defined above, and a kinetic
temperature that controls the velocity distribution of particles as
usual. Either constant volume or constant pressure algorithms can be
used.
The fix enforces a generalized ensemble in a single replica
only. Typically, this ideaology is combined with replica
exchange with replicas differing by {lambda} only
for simplicity, but this is not required. A multi-replica
simulation can be run within the LAMMPS environment using the
"grem"_temper_grem.html command. This utilizes LAMMPS partition
mode and requires the number of available processors be
on the order of the number of desired replicas. A 100-replica
simulation would require at least 100 processors (1 per world
at minimum). If a many replicas are needed on a small number of
processors, multi-replica runs can be run outside of LAMMPS.
An example of this can be found in examples/USER/misc/grem and
has no limit on the number of replicas per processor. However,
this is very inefficient and error prone and should be avoided
if possible.
In general, defining the generalized ensembles is unique for every
system. When starting a many-replica simulation without any knowledge of
the underlying microcanonical temperature, there are several tricks we
have utilized to optimize the process. Choosing a less-steep {eta}
yields broader distributions, requiring fewer replicas to map the
microcanonical temperature. While this likely struggles from the same
sampling problems gREM was built to avoid, it provides quick insight to
Ts. Initially using an evenly-spaced {lambda} distribution identifies
regions where small changes in enthalpy lead to large temperature
changes. Replicas are easily added where needed.
:line
[Restart, fix_modify, output, run start/stop, minimize info:]
No information about this fix is written to "binary restart
files"_restart.html.
The "thermo_modify"_thermo_modify.html {press} option is supported
by this fix to add the rescaled kinetic pressure as part of
"thermodynamic output"_thermo_style.html.
[Restrictions:]
This fix is part of the USER-MISC 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.
[Related commands:]
"grem"_temper_grem.html, "fix nvt"_fix_nh.html, "fix npt"_fix_nh.html, "thermo_modify"_thermo_modify.html
[Default:] none
:line
:link(Kim)
[(Kim)] Kim, Keyes, Straub, J. Chem. Phys., 132, 224107 (2010).
:link(Malolepsza)
[(Malolepsza)] Malolepsza, Secor, Keyes, J. Phys. Chem. B. 119 (42),
13379-13384 (2015).

View File

@ -48,6 +48,7 @@ Fixes :h1
fix_gld
fix_gle
fix_gravity
fix_grem
fix_halt
fix_heat
fix_imd

View File

@ -172,6 +172,7 @@ fix_gcmc.html
fix_gld.html
fix_gle.html
fix_gravity.html
fix_grem.html
fix_halt.html
fix_heat.html
fix_imd.html

101
doc/src/temper_grem.txt Normal file
View File

@ -0,0 +1,101 @@
"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
grem command :h3
[Syntax:]
grem N M lambda fix-ID thermostat-ID seed1 seed2 index :pre
N = total # of timesteps to run
M = attempt a tempering swap every this many steps
lambda = initial lambda for this ensemble
fix-ID = ID of fix_grem
thermostat-ID = ID of the thermostat that controls kinetic temperature
seed1 = random # seed used to decide on adjacent temperature to partner with
seed2 = random # seed for Boltzmann factor in Metropolis swap
index = which temperature (0 to N-1) I am simulating (optional) :ul
[Examples:]
grem 100000 1000 ${lambda} fxgREM fxnvt 0 58728
grem 40000 100 ${lambda} fxgREM fxnpt 0 32285 ${walkers} :pre
[Description:]
Run a parallel tempering or replica exchange simulation in LAMMPS partition
mode using multiple generalized replicas (ensembles) of a system defined by
"fix grem"_fix_grem.html. Two or more replicas must be used.
This command is a modification of the "temper"_temper.html command and has
the same dependencies, restraints, and input variables which are discussed
there in greater detail.
Instead of temperature, this command performs replica exchanges in lambda
as per the generalized ensemble enforced by "fix grem"_fix_grem.html.
The desired lambda is specified by {lambda}, which is typically a
variable previously set in the input script, so that each partition is
assigned a different temperature. See the "variable"_variable.html
command for more details. For example:
variable lambda world 400 420 440 460
fix fxnvt all nvt temp 300.0 300.0 100.0
fix fxgREM all grem ${lambda} -0.05 -50000 fxnvt
temper 100000 100 ${lambda} fxgREM fxnvt 3847 58382 :pre
would define 4 lambdas with constant kinetic temperature but unique
generalized temperature, and assign one of them to
"fix grem"_fix_grem.html used by each replica, and to the grem command.
As the gREM simulation runs for {N} timesteps, a swap
between adjacent ensembles will be attempted every {M} timesteps. If
{seed1} is 0, then the swap attempts will alternate between odd and
even pairings. If {seed1} is non-zero then it is used as a seed in a
random number generator to randomly choose an odd or even pairing each
time. Each attempted swap of temperatures is either accepted or
rejected based on a Metropolis criterion, derived for gREM by
"(Kim)"_#Kim, which uses {seed2} in the random number generator.
File management works identical to the "temper"_temper.html command.
Dump files created by this fix contain continuous trajectories and
require post-processing to obtain per-replica information.
The last argument {index} in the grem command is optional and is
used when restarting a run from a set of restart files (one
for each replica) which had previously swapped to new lambda.
This is done using a variable. For example if the
log file listed the following for a simulation with 5 replicas:
500000 2 4 0 1 3 :pre
then a setting of
variable walkers world 2 4 0 1 3 :pre
would be used to restart the run with a grem command like the
example above with ${walkers} as the last argument. This functionality
is identical to "temper"_temper.html.
:line
[Restrictions:]
This command can only be used if LAMMPS was built with the USER-MISC
package. See the "Making LAMMPS"_Section_start.html#start_3 section
for more info on packages.
This command must be used with "fix grem"_fix_grem.html.
[Related commands:]
"fix grem"_fix_grem.html, "temper"_temper.html, "variable"_variable.html
[Default:] none
:link(Kim)
[(Kim)] Kim, Keyes, Straub, J. Chem. Phys., 132, 224107 (2010).

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

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,15 @@
#/bin/bash
for i in $(ls -d [0-9]*)
do
rm -f $i/final*
rm -f $i/log*
rm -f $i/ent*
rm -f $i/output
cp $i/restart.init $i/restart_file
done
echo 1 > lastexchange
cp walker.bkp lastwalker
exit 0

View File

@ -0,0 +1,169 @@
#!/usr/bin/env python2.7
import os, sys
from numpy import *
import numpy.random
### Runs replica exchange with gREM (fix grem) for unlimited number of replicas on a set number of processors. This script is inefficient, but necessary if wanting to run with hundreds of replicas on relatively few number of procs.
### read number of processors from the command line
nproc = int(sys.argv[1])
### path to simulation directory
path = os.getcwd()
### path to LAMMPS executable
lmp = sys.argv[2]
### LAMMPS input name
inp = sys.argv[3]
### define pressure for simulations (0 if const V)
pressure = 0
### some constants for gREM, must match with LAMMPS input file!
H = -30000
eta = -0.01
#kB = 0.000086173324 # eV (metal)
kB = 0.0019872 # kcal/mol (real)
### define lambdas - script assumes that there are already existing directories with all files necessary to run
lambdas=[400,405,410,415,420,425]
ll = len(lambdas)
### define number of exchanges
starting_ex = int(loadtxt("lastexchange"))
how_many_ex = 5
max_exchange = starting_ex+how_many_ex
### array with walkers
walker = loadtxt("lastwalker")
### initiate array with enthalpies
enthalpy = zeros(ll)
aver_enthalpy = zeros(ll)
for exchange in arange(starting_ex,max_exchange):
print "run", exchange
for l in range(ll):
#print "replica", l
os.chdir(path+"/%s" % lambdas[l])
#os.system("cp restart_file restart_file%d" % exchange)
if (nproc > 1):
os.system("mpirun -np %d " + lmp + " -in ../" + inp + " -var lambda %g -var eta %g -var enthalpy %g > output" % (nproc, lambdas[l], eta, H))
if (nproc == 1):
os.system(lmp + " -in ../" + inp + " -var lambda %g -var eta %g -var enthalpy %g > output" % (lambdas[l], eta, H))
os.system("grep -v '[a-zA-Z]' output | awk '{if(NF==6 && NR>19)print $0}' | awk '{print $3}' >ent")
enthalpy[l] = os.popen("tail -n 1 ent").read()
ee = loadtxt("ent")
aver_enthalpy[l] = mean(ee[-1])
# os.system("mv dump.dcd dump%d.dcd" % exchange)
os.system("mv log.lammps log%d.lammps" % exchange)
# os.system("rm output")
os.system("mv final_restart_file final_restart_file%d" % exchange)
os.system("mv ent ent%d" % exchange)
os.system("bzip2 log%d.lammps ent%d" % (exchange,exchange))
os.system("cp final_restart_file%d restart_file" % exchange)
### replicas will be exchanged based on enthalpy order, not replicas order (termostat order)
#entalpy_sorted_indices = enthalpy.argsort()
aver_entalpy_sorted_indices = aver_enthalpy.argsort()
### choose pair of replicas for exchange attempt based on enthalpy order
pp = random.random_integers(0,ll-2)
first = aver_entalpy_sorted_indices[pp]
second = aver_entalpy_sorted_indices[pp+1]
#if (first>second):
# tmp = first
# first = second
# second = tmp
print "pair1:", first, second
### calculate weights for exchange criterion
w1 = log(lambdas[first]+eta*(enthalpy[first]-1*H))
w2 = log(lambdas[first]+eta*(enthalpy[second]-1*H))
w3 = log(lambdas[second]+eta*(enthalpy[first]-1*H))
w4 = log(lambdas[second]+eta*(enthalpy[second]-1*H))
weight = (w4-w3+w1-w2)/eta/kB
### generate randon number for exchange criterion and calc its log
LOGRANDNUM = log(random.random())
### wyzeruj warunki
compare1 = 0
compare2 = 0
if (weight>0):
compare1 = 1
if (weight>LOGRANDNUM):
compare2 = 1
### exchange restart files if exchange condition is satisfied
if (compare1>0 or compare2>0):
print "exchange1 accepted for pair", first, second, lambdas[first], lambdas[second], "with compares as", compare1, compare2, "weight as", weight, "and lograndnum", LOGRANDNUM
os.system("cp %s/%s/final_restart_file%d %s/%s/restart_file" % (path,lambdas[first],exchange,path,lambdas[second]))
os.system("cp %s/%s/final_restart_file%d %s/%s/restart_file" % (path,lambdas[second],exchange,path,lambdas[first]))
### update walkers
tmp1=walker[first]
tmp2=walker[second]
walker[first]=tmp2
walker[second]=tmp1
else:
print "exchange1 not accepted for pair", first, second, lambdas[first], lambdas[second], "with compares as", compare1, compare2, "weight as", weight, "and lograndnum", LOGRANDNUM
### choose again pair of replicas for exchange attempt based on enthalpy order
### but make sure this pair is different than the first pair
if_different = 0
while if_different<1:
pp2 = random.random_integers(0,ll-2)
third = aver_entalpy_sorted_indices[pp2]
fourth = aver_entalpy_sorted_indices[pp2+1]
if (third!=first and third!=second and third!=aver_entalpy_sorted_indices[pp-1]):
if_different = 1
print "pair2:", third, fourth
### calculate weights for exchange criterion
w1 = log(lambdas[third]+eta*(enthalpy[third]-1*H))
w2 = log(lambdas[third]+eta*(enthalpy[fourth]-1*H))
w3 = log(lambdas[fourth]+eta*(enthalpy[third]-1*H))
w4 = log(lambdas[fourth]+eta*(enthalpy[fourth]-1*H))
weight = (w4-w3+w1-w2)/eta/kB
### generate randon number for exchange criterion and calc its log
LOGRANDNUM = log(random.random())
### wyzeruj warunki
compare1 = 0
compare2 = 0
if (weight>0):
compare1 = 1
if (weight>LOGRANDNUM):
compare2 = 1
### exchange restart files if exchange condition is satisfied
if (compare1>0 or compare2>0):
print "exchange2 accepted for pair", third, fourth, lambdas[third], lambdas[fourth], "with compares as", compare1, compare2, "weight as", weight, "and lograndnum", LOGRANDNUM
os.system("cp %s/%s/final_restart_file%d %s/%s/restart_file" % (path,lambdas[third],exchange,path,lambdas[fourth]))
os.system("cp %s/%s/final_restart_file%d %s/%s/restart_file" % (path,lambdas[fourth],exchange,path,lambdas[third]))
### update walkers
tmp1=walker[third]
tmp2=walker[fourth]
walker[third]=tmp2
walker[fourth]=tmp1
else:
print "exchange2 not accepted for pair", third, fourth, lambdas[third], lambdas[fourth], "with compares as", compare1, compare2, "weight as", weight, "and lograndnum", LOGRANDNUM
#print "walkers:", walker
print "".join(["%d " % x for x in walker])
sys.stdout.flush()
lastwalker = open(path + "/lastwalker", "w")
lastwalker.write("".join(["%d " % w for w in walker]))
lastwalker.close()
lastexchange = open(path + "/lastexchange", "w")
lastexchange.write("%d" % (exchange+1))
lastexchange.close()

View File

@ -0,0 +1,25 @@
# LJ particles
variable T0 index 300.0
variable press index 0.0
variable lambda index 400.0
variable eta index -0.01
variable enthalpy index -30000.0
units real
atom_style full
pair_style lj/cut 5.0
read_data "restart_file"
thermo 10
thermo_style custom step temp pe etotal press vol
velocity all create ${T0} 12427
timestep 1.0
fix fxnvt all npt temp ${T0} ${T0} 1000.0 iso ${press} ${press} 10000.0
fix fxgREM all grem ${lambda} ${eta} ${enthalpy} fxnvt
thermo_modify press fxgREM_press
run 10000
write_data final_restart_file

View File

@ -0,0 +1,12 @@
#!/bin/sh
NPROCS=1
if [ $# -gt 0 ]; then
NPROCS=$1
fi
bash ./clean.sh
python ./double-re-short.py $NPROCS $HOME/compile/lammps-icms/src/lmp_omp in.gREM > total_output.$NPROCS
exit 0

View File

@ -0,0 +1 @@
0 1 2 3 4 5

View File

@ -0,0 +1,21 @@
# LJ particles
variable T0 equal 300.0
variable press equal 0.0
units real
atom_style full
pair_style lj/cut 5.0
read_data "lj.data"
thermo 10
thermo_style custom step temp pe etotal press vol
timestep 1.0
fix fxnpt all npt temp ${T0} ${T0} 1000.0 iso ${press} ${press} 10000.0
fix fxgREM all grem 400 -.01 -30000 fxnpt
thermo_modify press fxgREM_press
run 1000
#write_data lj-out.data

View File

@ -0,0 +1,20 @@
# LJ particles
variable T0 equal 300.0
variable press equal 0.0
units real
atom_style full
pair_style lj/cut 5.0
read_data "lj.data"
thermo 10
thermo_style custom step temp pe etotal press vol
timestep 1.0
fix fxnvt all nvt temp ${T0} ${T0} 1000.0
fix fxgREM all grem 400 -.01 -30000 fxnvt
run 1000
#write_data lj-out.data

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,176 @@
LAMMPS (9 Nov 2016)
using 1 OpenMP thread(s) per MPI task
# LJ particles
variable T0 equal 300.0
variable press equal 0.0
units real
atom_style full
pair_style lj/cut 5.0
read_data "lj.data"
orthogonal box = (1.06874 1.06874 1.06874) to (23.9313 23.9313 23.9313)
1 by 1 by 1 MPI processor grid
reading atoms ...
500 atoms
reading velocities ...
500 velocities
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
thermo 10
thermo_style custom step temp pe etotal press vol
timestep 1.0
fix fxnvt all npt temp ${T0} ${T0} 1000.0 iso ${press} ${press} 10000.0
fix fxnvt all npt temp 300 ${T0} 1000.0 iso ${press} ${press} 10000.0
fix fxnvt all npt temp 300 300 1000.0 iso ${press} ${press} 10000.0
fix fxnvt all npt temp 300 300 1000.0 iso 0 ${press} 10000.0
fix fxnvt all npt temp 300 300 1000.0 iso 0 0 10000.0
fix fxgREM all grem 400 -.01 -30000 fxnvt
thermo_modify press fxgREM_press
run 1000
Neighbor list info ...
1 neighbor list requests
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 7
ghost atom cutoff = 7
binsize = 3.5 -> bins = 7 7 7
Memory usage per processor = 5.37943 Mbytes
Step Temp PotEng TotEng Press Volume
0 305.69499 -3177.6423 -2722.9442 -91.741776 11950.115
10 312.30124 -3182.2257 -2717.7013 -203.95075 11950.113
20 314.94567 -3186.456 -2717.9982 -265.56737 11950.108
30 312.229 -3183.7641 -2719.3472 -196.90499 11950.097
40 305.94068 -3180.7085 -2725.6449 -92.562221 11950.083
50 300.42281 -3176.5838 -2729.7277 10.896769 11950.066
60 299.16747 -3174.1939 -2729.205 50.094171 11950.05
70 301.65965 -3176.0918 -2727.396 0.096901939 11950.035
80 304.77876 -3178.2699 -2724.9346 -64.001022 11950.019
90 305.60598 -3178.9517 -2724.386 -93.672879 11950.003
100 303.8005 -3177.5156 -2725.6354 -74.516709 11949.985
110 300.86776 -3175.4773 -2727.9593 -34.22655 11949.965
120 298.70177 -3175.6488 -2731.3526 -19.014898 11949.944
130 298.39686 -3176.3792 -2732.5365 -21.293245 11949.923
140 300.00669 -3177.7032 -2731.466 -40.992937 11949.902
150 301.85665 -3178.1312 -2729.1423 -45.715505 11949.88
160 301.20597 -3177.3218 -2729.3007 -10.104082 11949.857
170 297.01134 -3172.7462 -2730.9643 99.298381 11949.833
180 291.279 -3168.3513 -2735.0958 219.47549 11949.812
190 287.13954 -3165.1287 -2738.0304 309.36947 11949.796
200 286.57735 -3165.2951 -2739.033 323.96954 11949.786
210 289.83941 -3167.8245 -2736.7103 271.77305 11949.783
220 296.12858 -3171.8054 -2731.3366 172.4056 11949.785
230 303.82424 -3176.3108 -2724.3952 56.711479 11949.791
240 309.95738 -3180.9789 -2719.9408 -40.992898 11949.798
250 312.0405 -3182.3473 -2718.2107 -57.591676 11949.805
260 309.65444 -3181.0587 -2720.4712 3.3540332 11949.81
270 304.40001 -3176.5798 -2723.8078 130.77028 11949.816
280 298.65985 -3174.1505 -2729.9166 237.63562 11949.825
290 294.78709 -3170.9701 -2732.4966 326.94924 11949.838
300 294.03216 -3169.9567 -2732.6062 349.85486 11949.859
310 296.44397 -3172.8519 -2731.914 284.80897 11949.886
320 301.41027 -3175.9697 -2727.6447 179.4647 11949.92
330 307.88911 -3181.2615 -2723.2998 24.702414 11949.957
340 314.73138 -3186.0047 -2717.8656 -132.6263 11949.995
350 320.55591 -3187.8509 -2711.0483 -245.88468 11950.031
360 323.50274 -3188.9994 -2707.8136 -314.73676 11950.062
370 321.61539 -3187.1233 -2708.7448 -293.17446 11950.086
380 314.37275 -3181.484 -2713.8784 -169.00448 11950.104
390 303.54884 -3174.1675 -2722.6616 12.923999 11950.119
400 293.40432 -3167.0348 -2730.6181 187.6624 11950.135
410 288.46351 -3165.273 -2736.2054 252.20051 11950.154
420 290.31387 -3168.604 -2736.7841 193.73816 11950.178
430 296.35519 -3173.09 -2732.2841 81.521847 11950.207
440 301.92973 -3175.4344 -2726.3368 -1.8329439 11950.237
450 303.76205 -3176.777 -2724.9539 -35.002096 11950.267
460 301.71619 -3174.2731 -2725.4932 14.977875 11950.296
470 298.92404 -3172.9921 -2728.3652 64.224747 11950.326
480 298.80164 -3172.5329 -2728.0881 82.781347 11950.358
490 302.71589 -3175.3703 -2725.1034 27.223049 11950.39
500 309.10665 -3179.3013 -2719.5285 -65.460658 11950.424
510 314.36408 -3183.2854 -2715.6927 -151.19245 11950.456
520 315.71154 -3183.5328 -2713.9358 -163.19151 11950.485
530 313.31886 -3182.2521 -2716.214 -125.5741 11950.511
540 309.81847 -3178.9358 -2718.1043 -55.55841 11950.534
550 308.29687 -3177.837 -2719.2688 -24.39371 11950.556
560 308.75927 -3176.3265 -2717.0705 0.93689833 11950.578
570 307.52811 -3175.8145 -2718.3897 35.502429 11950.6
580 301.75074 -3173.1208 -2724.2894 136.29625 11950.622
590 292.37743 -3165.5806 -2730.6913 319.75957 11950.648
600 283.57627 -3159.8617 -2738.0635 471.28045 11950.68
610 279.85172 -3157.4557 -2741.1975 530.72699 11950.722
620 283.40879 -3160.5911 -2739.042 455.28104 11950.775
630 292.53718 -3166.3125 -2731.1856 296.63465 11950.838
640 302.81112 -3173.3096 -2722.901 113.80844 11950.907
650 309.83321 -3179.3684 -2718.515 -26.499431 11950.978
660 312.1283 -3182.7335 -2718.4663 -89.363745 11951.049
670 311.16363 -3181.867 -2719.0347 -69.370989 11951.118
680 308.51041 -3180.6869 -2721.801 -25.972987 11951.186
690 304.64393 -3176.8751 -2723.7403 56.592367 11951.254
700 300.24456 -3175.4797 -2728.8887 112.34442 11951.323
710 296.35785 -3172.9705 -2732.1607 168.18009 11951.394
720 293.78145 -3172.1065 -2735.1289 182.81082 11951.468
730 293.25707 -3170.8715 -2734.6738 171.04236 11951.547
740 295.33219 -3172.9109 -2733.6266 91.351362 11951.629
750 299.69136 -3175.2574 -2729.4892 -16.266404 11951.713
760 305.2281 -3177.9836 -2723.9799 -137.30615 11951.796
770 310.59309 -3182.7053 -2720.7216 -272.72961 11951.877
780 314.65573 -3183.4212 -2715.3947 -341.231 11951.952
790 316.48606 -3185.44 -2714.691 -388.53602 11952.02
800 315.15897 -3186.846 -2718.0709 -384.28316 11952.08
810 310.43559 -3183.6648 -2721.9154 -282.61999 11952.133
820 303.22265 -3178.464 -2727.4433 -121.47565 11952.179
830 295.36843 -3175.4771 -2736.1389 33.066504 11952.223
840 288.69698 -3169.5813 -2740.1664 216.10697 11952.268
850 283.82649 -3165.7822 -2743.6118 359.56896 11952.317
860 280.04102 -3162.8228 -2746.283 475.61942 11952.374
870 277.10059 -3159.6212 -2747.4551 572.5432 11952.441
880 275.76549 -3158.2545 -2748.0743 616.43304 11952.52
890 276.82327 -3158.9703 -2747.2166 596.08147 11952.612
900 280.72135 -3162.0637 -2744.5119 506.33695 11952.716
910 287.1035 -3167.4388 -2740.3941 356.68688 11952.831
920 294.28041 -3171.6218 -2733.902 206.06394 11952.953
930 300.36009 -3173.9046 -2727.1418 88.047911 11953.08
940 303.86761 -3175.5599 -2723.5798 7.6846808 11953.209
950 304.42957 -3176.0831 -2723.2672 -25.15496 11953.339
960 303.13982 -3176.0534 -2725.1559 -28.715178 11953.467
970 302.30166 -3176.9758 -2727.325 -43.264668 11953.596
980 303.93331 -3178.9891 -2726.9114 -88.434034 11953.723
990 307.36223 -3180.7316 -2723.5535 -145.46208 11953.849
1000 310.09574 -3181.101 -2719.8571 -180.39125 11953.972
Loop time of 0.307225 on 1 procs for 1000 steps with 500 atoms
Performance: 281.227 ns/day, 0.085 hours/ns, 3254.944 timesteps/s
99.6% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.25351 | 0.25351 | 0.25351 | 0.0 | 82.52
Bond | 7.1526e-05 | 7.1526e-05 | 7.1526e-05 | 0.0 | 0.02
Neigh | 0.0042093 | 0.0042093 | 0.0042093 | 0.0 | 1.37
Comm | 0.010211 | 0.010211 | 0.010211 | 0.0 | 3.32
Output | 0.0013611 | 0.0013611 | 0.0013611 | 0.0 | 0.44
Modify | 0.033891 | 0.033891 | 0.033891 | 0.0 | 11.03
Other | | 0.003969 | | | 1.29
Nlocal: 500 ave 500 max 500 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1610 ave 1610 max 1610 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 14765 ave 14765 max 14765 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 14765
Ave neighs/atom = 29.53
Ave special neighs/atom = 0
Neighbor list builds = 5
Dangerous builds = 0
#write_data lj-out.data
Total wall time: 0:00:00

View File

@ -0,0 +1,176 @@
LAMMPS (9 Nov 2016)
using 1 OpenMP thread(s) per MPI task
# LJ particles
variable T0 equal 300.0
variable press equal 0.0
units real
atom_style full
pair_style lj/cut 5.0
read_data "lj.data"
orthogonal box = (1.06874 1.06874 1.06874) to (23.9313 23.9313 23.9313)
1 by 2 by 2 MPI processor grid
reading atoms ...
500 atoms
reading velocities ...
500 velocities
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
thermo 10
thermo_style custom step temp pe etotal press vol
timestep 1.0
fix fxnvt all npt temp ${T0} ${T0} 1000.0 iso ${press} ${press} 10000.0
fix fxnvt all npt temp 300 ${T0} 1000.0 iso ${press} ${press} 10000.0
fix fxnvt all npt temp 300 300 1000.0 iso ${press} ${press} 10000.0
fix fxnvt all npt temp 300 300 1000.0 iso 0 ${press} 10000.0
fix fxnvt all npt temp 300 300 1000.0 iso 0 0 10000.0
fix fxgREM all grem 400 -.01 -30000 fxnvt
thermo_modify press fxgREM_press
run 1000
Neighbor list info ...
1 neighbor list requests
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 7
ghost atom cutoff = 7
binsize = 3.5 -> bins = 7 7 7
Memory usage per processor = 5.34276 Mbytes
Step Temp PotEng TotEng Press Volume
0 305.69499 -3177.6423 -2722.9442 -91.741776 11950.115
10 312.30124 -3182.2257 -2717.7013 -203.95075 11950.113
20 314.94567 -3186.456 -2717.9982 -265.56737 11950.108
30 312.229 -3183.7641 -2719.3472 -196.90499 11950.097
40 305.94068 -3180.7085 -2725.6449 -92.562221 11950.083
50 300.42281 -3176.5838 -2729.7277 10.896769 11950.066
60 299.16747 -3174.1939 -2729.205 50.094171 11950.05
70 301.65965 -3176.0918 -2727.396 0.096901939 11950.035
80 304.77876 -3178.2699 -2724.9346 -64.001022 11950.019
90 305.60598 -3178.9517 -2724.386 -93.672879 11950.003
100 303.8005 -3177.5156 -2725.6354 -74.516709 11949.985
110 300.86776 -3175.4773 -2727.9593 -34.22655 11949.965
120 298.70177 -3175.6488 -2731.3526 -19.014898 11949.944
130 298.39686 -3176.3792 -2732.5365 -21.293245 11949.923
140 300.00669 -3177.7032 -2731.466 -40.992937 11949.902
150 301.85665 -3178.1312 -2729.1423 -45.715505 11949.88
160 301.20597 -3177.3218 -2729.3007 -10.104082 11949.857
170 297.01134 -3172.7462 -2730.9643 99.298381 11949.833
180 291.279 -3168.3513 -2735.0958 219.47549 11949.812
190 287.13954 -3165.1287 -2738.0304 309.36947 11949.796
200 286.57735 -3165.2951 -2739.033 323.96954 11949.786
210 289.83941 -3167.8245 -2736.7103 271.77305 11949.783
220 296.12858 -3171.8054 -2731.3366 172.4056 11949.785
230 303.82424 -3176.3108 -2724.3952 56.711479 11949.791
240 309.95738 -3180.9789 -2719.9408 -40.992898 11949.798
250 312.0405 -3182.3473 -2718.2107 -57.591676 11949.805
260 309.65444 -3181.0587 -2720.4712 3.3540332 11949.81
270 304.40001 -3176.5798 -2723.8078 130.77028 11949.816
280 298.65985 -3174.1505 -2729.9166 237.63562 11949.825
290 294.78709 -3170.9701 -2732.4966 326.94924 11949.838
300 294.03216 -3169.9567 -2732.6062 349.85486 11949.859
310 296.44397 -3172.8519 -2731.914 284.80897 11949.886
320 301.41027 -3175.9697 -2727.6447 179.4647 11949.92
330 307.88911 -3181.2615 -2723.2998 24.702414 11949.957
340 314.73138 -3186.0047 -2717.8656 -132.6263 11949.995
350 320.55591 -3187.8509 -2711.0483 -245.88468 11950.031
360 323.50274 -3188.9994 -2707.8136 -314.73676 11950.062
370 321.61539 -3187.1233 -2708.7448 -293.17446 11950.086
380 314.37275 -3181.484 -2713.8784 -169.00448 11950.104
390 303.54884 -3174.1675 -2722.6616 12.923999 11950.119
400 293.40432 -3167.0348 -2730.6181 187.6624 11950.135
410 288.46351 -3165.273 -2736.2054 252.20051 11950.154
420 290.31387 -3168.604 -2736.7841 193.73816 11950.178
430 296.35519 -3173.09 -2732.2841 81.521847 11950.207
440 301.92973 -3175.4344 -2726.3368 -1.8329439 11950.237
450 303.76205 -3176.777 -2724.9539 -35.002096 11950.267
460 301.71619 -3174.2731 -2725.4932 14.977875 11950.296
470 298.92404 -3172.9921 -2728.3652 64.224747 11950.326
480 298.80164 -3172.5329 -2728.0881 82.781347 11950.358
490 302.71589 -3175.3703 -2725.1034 27.223049 11950.39
500 309.10665 -3179.3013 -2719.5285 -65.460658 11950.424
510 314.36408 -3183.2854 -2715.6927 -151.19245 11950.456
520 315.71154 -3183.5328 -2713.9358 -163.19151 11950.485
530 313.31886 -3182.2521 -2716.214 -125.5741 11950.511
540 309.81847 -3178.9358 -2718.1043 -55.55841 11950.534
550 308.29687 -3177.837 -2719.2688 -24.39371 11950.556
560 308.75927 -3176.3265 -2717.0705 0.93689833 11950.578
570 307.52811 -3175.8145 -2718.3897 35.502429 11950.6
580 301.75074 -3173.1208 -2724.2894 136.29625 11950.622
590 292.37743 -3165.5806 -2730.6913 319.75957 11950.648
600 283.57627 -3159.8617 -2738.0635 471.28045 11950.68
610 279.85172 -3157.4557 -2741.1975 530.72699 11950.722
620 283.40879 -3160.5911 -2739.042 455.28104 11950.775
630 292.53718 -3166.3125 -2731.1856 296.63465 11950.838
640 302.81112 -3173.3096 -2722.901 113.80844 11950.907
650 309.83321 -3179.3684 -2718.515 -26.499431 11950.978
660 312.1283 -3182.7335 -2718.4663 -89.363745 11951.049
670 311.16363 -3181.867 -2719.0347 -69.370989 11951.118
680 308.51041 -3180.6869 -2721.801 -25.972987 11951.186
690 304.64393 -3176.8751 -2723.7403 56.592367 11951.254
700 300.24456 -3175.4797 -2728.8887 112.34442 11951.323
710 296.35785 -3172.9705 -2732.1607 168.18009 11951.394
720 293.78145 -3172.1065 -2735.1289 182.81082 11951.468
730 293.25707 -3170.8715 -2734.6738 171.04236 11951.547
740 295.33219 -3172.9109 -2733.6266 91.351362 11951.629
750 299.69136 -3175.2574 -2729.4892 -16.266404 11951.713
760 305.2281 -3177.9836 -2723.9799 -137.30615 11951.796
770 310.59309 -3182.7053 -2720.7216 -272.72961 11951.877
780 314.65573 -3183.4212 -2715.3947 -341.231 11951.952
790 316.48606 -3185.44 -2714.691 -388.53602 11952.02
800 315.15897 -3186.846 -2718.0709 -384.28316 11952.08
810 310.43559 -3183.6648 -2721.9154 -282.61999 11952.133
820 303.22265 -3178.464 -2727.4433 -121.47565 11952.179
830 295.36843 -3175.4771 -2736.1389 33.066504 11952.223
840 288.69698 -3169.5813 -2740.1664 216.10697 11952.268
850 283.82649 -3165.7822 -2743.6118 359.56896 11952.317
860 280.04102 -3162.8228 -2746.283 475.61942 11952.374
870 277.10059 -3159.6212 -2747.4551 572.5432 11952.441
880 275.76549 -3158.2545 -2748.0743 616.43304 11952.52
890 276.82327 -3158.9703 -2747.2166 596.08147 11952.612
900 280.72135 -3162.0637 -2744.5119 506.33695 11952.716
910 287.1035 -3167.4388 -2740.3941 356.68688 11952.831
920 294.28041 -3171.6218 -2733.902 206.06394 11952.953
930 300.36009 -3173.9046 -2727.1418 88.047911 11953.08
940 303.86761 -3175.5599 -2723.5798 7.6846808 11953.209
950 304.42957 -3176.0831 -2723.2672 -25.15496 11953.339
960 303.13982 -3176.0534 -2725.1559 -28.715178 11953.467
970 302.30166 -3176.9758 -2727.325 -43.264668 11953.596
980 303.93331 -3178.9891 -2726.9114 -88.434034 11953.723
990 307.36223 -3180.7316 -2723.5535 -145.46208 11953.849
1000 310.09574 -3181.101 -2719.8571 -180.39125 11953.972
Loop time of 0.154208 on 4 procs for 1000 steps with 500 atoms
Performance: 560.281 ns/day, 0.043 hours/ns, 6484.730 timesteps/s
98.1% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.072079 | 0.074846 | 0.079666 | 1.1 | 48.54
Bond | 5.7936e-05 | 6.634e-05 | 8.1062e-05 | 0.1 | 0.04
Neigh | 0.0010812 | 0.0012064 | 0.0012748 | 0.2 | 0.78
Comm | 0.032452 | 0.037544 | 0.04076 | 1.6 | 24.35
Output | 0.0018461 | 0.0020589 | 0.0026393 | 0.7 | 1.34
Modify | 0.032085 | 0.032688 | 0.033361 | 0.3 | 21.20
Other | | 0.005799 | | | 3.76
Nlocal: 125 ave 127 max 123 min
Histogram: 1 0 1 0 0 0 0 1 0 1
Nghost: 870.5 ave 882 max 862 min
Histogram: 1 1 0 0 0 0 1 0 0 1
Neighs: 3691.25 ave 3807 max 3563 min
Histogram: 1 0 0 0 1 0 1 0 0 1
Total # of neighbors = 14765
Ave neighs/atom = 29.53
Ave special neighs/atom = 0
Neighbor list builds = 5
Dangerous builds = 0
#write_data lj-out.data
Total wall time: 0:00:00

View File

@ -0,0 +1,173 @@
LAMMPS (9 Nov 2016)
using 1 OpenMP thread(s) per MPI task
# LJ particles
variable T0 equal 300.0
variable press equal 0.0
units real
atom_style full
pair_style lj/cut 5.0
read_data "lj.data"
orthogonal box = (1.06874 1.06874 1.06874) to (23.9313 23.9313 23.9313)
1 by 1 by 1 MPI processor grid
reading atoms ...
500 atoms
reading velocities ...
500 velocities
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
thermo 10
thermo_style custom step temp pe etotal press vol
timestep 1.0
fix fxnvt all nvt temp ${T0} ${T0} 1000.0
fix fxnvt all nvt temp 300 ${T0} 1000.0
fix fxnvt all nvt temp 300 300 1000.0
fix fxgREM all grem 400 -.01 -30000 fxnvt
run 1000
Neighbor list info ...
1 neighbor list requests
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 7
ghost atom cutoff = 7
binsize = 3.5 -> bins = 7 7 7
Memory usage per processor = 5.37943 Mbytes
Step Temp PotEng TotEng Press Volume
0 305.69499 -3177.6423 -2722.9442 883.58369 11950.115
10 312.30121 -3182.2257 -2717.7013 793.47811 11950.115
20 314.94553 -3186.4559 -2717.9983 738.74091 11950.115
30 312.22861 -3183.7638 -2719.3474 797.47978 11950.115
40 305.93987 -3180.7079 -2725.6455 881.30806 11950.115
50 300.42132 -3176.5828 -2729.7288 967.92042 11950.115
60 299.16487 -3174.1921 -2729.2071 1004.247 11950.115
70 301.65565 -3176.0891 -2727.3992 962.58134 11950.115
80 304.77334 -3178.2663 -2724.939 907.8946 11950.115
90 305.59929 -3178.9472 -2724.3914 879.91629 11950.115
100 303.79263 -3177.5103 -2725.6418 892.67631 11950.115
110 300.85863 -3175.4711 -2727.9667 923.44924 11950.115
120 298.69083 -3175.6415 -2731.3615 931.87518 11950.115
130 298.38415 -3176.3706 -2732.5468 928.88286 11950.115
140 299.99129 -3177.6935 -2731.4792 914.36783 11950.115
150 301.83869 -3178.121 -2729.1588 915.01407 11950.115
160 301.18834 -3177.3117 -2729.3169 947.45228 11950.115
170 296.99406 -3172.7363 -2730.9801 1042.6928 11950.115
180 291.25952 -3168.3407 -2735.1142 1144.5436 11950.115
190 287.1178 -3164.9847 -2737.9187 1223.4003 11950.115
200 286.552 -3165.2799 -2739.0555 1235.6703 11950.115
210 289.81033 -3167.8062 -2736.7353 1194.6672 11950.115
220 296.09616 -3171.7847 -2731.3641 1115.8799 11950.115
230 303.79176 -3176.2893 -2724.4221 1024.6471 11950.115
240 309.9273 -3180.9591 -2719.9657 945.55045 11950.115
250 312.0159 -3182.3307 -2718.2306 934.36956 11950.115
260 309.63264 -3181.0452 -2720.4901 986.77385 11950.115
270 304.38172 -3176.568 -2723.8233 1097.264 11950.115
280 298.64188 -3174.1384 -2729.9313 1186.2239 11950.115
290 294.76686 -3170.9562 -2732.5128 1264.247 11950.115
300 294.00805 -3169.8091 -2732.4944 1287.4001 11950.115
310 296.41801 -3172.834 -2731.9347 1229.5624 11950.115
320 301.38477 -3175.9514 -2727.6644 1140.8664 11950.115
330 307.86584 -3181.2442 -2723.3171 1007.1545 11950.115
340 314.7103 -3185.9891 -2717.8814 871.74528 11950.115
350 320.53954 -3187.8385 -2711.0602 776.85994 11950.115
360 323.49505 -3188.9927 -2707.8184 716.58062 11950.115
370 321.62077 -3187.1246 -2708.7381 731.01909 11950.115
380 314.39049 -3181.4931 -2713.8611 831.21057 11950.115
390 303.57079 -3174.1804 -2722.6419 978.62645 11950.115
400 293.42165 -3167.0452 -2730.6027 1122.3558 11950.115
410 288.46838 -3165.4071 -2736.3322 1171.8087 11950.115
420 290.30766 -3168.5988 -2736.7882 1122.5413 11950.115
430 296.34338 -3173.0824 -2732.2941 1030.2769 11950.115
440 301.92394 -3175.4307 -2726.3417 964.25387 11950.115
450 303.76745 -3176.9122 -2725.0811 934.49176 11950.115
460 301.72985 -3174.2821 -2725.4818 979.07605 11950.115
470 298.93736 -3173.0014 -2728.3548 1020.0482 11950.115
480 298.80912 -3172.803 -2728.3471 1036.6531 11950.115
490 302.72217 -3175.3764 -2725.1001 997.71146 11950.115
500 309.11393 -3179.3088 -2719.5253 925.81108 11950.115
510 314.37612 -3183.2961 -2715.6855 856.23748 11950.115
520 315.72767 -3183.547 -2713.926 847.70543 11950.115
530 313.34173 -3182.2695 -2716.1974 877.30842 11950.115
540 309.84312 -3178.9553 -2718.0871 936.69244 11950.115
550 308.3251 -3177.8582 -2719.248 963.93032 11950.115
560 308.79192 -3176.4834 -2717.1788 989.67643 11950.115
570 307.57194 -3175.8464 -2718.3565 1021.0494 11950.115
580 301.8035 -3173.1582 -2724.2483 1102.4893 11950.115
590 292.43425 -3165.751 -2730.7772 1254.7815 11950.115
600 283.62905 -3159.8987 -2738.022 1381.0608 11950.115
610 279.90122 -3157.49 -2741.1581 1431.0028 11950.115
620 283.4582 -3160.756 -2739.1334 1367.7385 11950.115
630 292.58866 -3166.3469 -2731.1435 1241.1194 11950.115
640 302.86585 -3173.4778 -2722.9878 1089.7342 11950.115
650 309.89252 -3179.4078 -2718.4662 972.6359 11950.115
660 312.19165 -3182.7754 -2718.414 916.62037 11950.115
670 311.2287 -3181.9102 -2718.9811 933.79804 11950.115
680 308.57852 -3180.7312 -2721.7441 969.24936 11950.115
690 304.71609 -3176.9196 -2723.6775 1040.2699 11950.115
700 300.31995 -3175.5245 -2728.8213 1082.845 11950.115
710 296.43537 -3173.0166 -2732.0915 1127.4487 11950.115
720 293.86692 -3172.1582 -2735.0535 1135.0215 11950.115
730 293.35611 -3170.9335 -2734.5885 1122.9143 11950.115
740 295.44861 -3172.9862 -2733.5288 1050.995 11950.115
750 299.82732 -3175.3467 -2729.3763 958.31462 11950.115
760 305.37987 -3178.216 -2723.9866 854.1946 11950.115
770 310.75394 -3182.8127 -2720.5898 737.72668 11950.115
780 314.81395 -3183.7905 -2715.5286 679.74198 11950.115
790 316.63339 -3185.8028 -2714.8346 638.48871 11950.115
800 315.2894 -3186.9345 -2717.9654 641.53256 11950.115
810 310.54289 -3183.7383 -2721.8293 728.51241 11950.115
820 303.31439 -3178.7897 -2727.6326 864.45674 11950.115
830 295.46125 -3175.5387 -2736.0625 997.72969 11950.115
840 288.802 -3169.6502 -2740.0791 1160.6622 11950.115
850 283.94785 -3165.8605 -2743.5096 1289.55 11950.115
860 280.17501 -3163.0381 -2746.299 1392.8854 11950.115
870 277.2456 -3159.8429 -2747.4611 1481.3899 11950.115
880 275.93123 -3158.3584 -2747.9316 1523.5374 11950.115
890 277.0215 -3159.2285 -2747.18 1506.1558 11950.115
900 280.96237 -3162.483 -2744.5728 1428.4183 11950.115
910 287.37962 -3167.6183 -2740.1628 1303.0268 11950.115
920 294.56731 -3171.6765 -2733.5299 1177.748 11950.115
930 300.63273 -3174.0842 -2726.9158 1078.7393 11950.115
940 304.10943 -3175.9847 -2723.645 1007.7154 11950.115
950 304.64845 -3176.6263 -2723.4848 976.37917 11950.115
960 303.36343 -3176.4694 -2725.2393 971.40749 11950.115
970 302.57138 -3177.5541 -2727.5021 954.01115 11950.115
980 304.2593 -3179.2101 -2726.6475 919.74949 11950.115
990 307.69959 -3180.9631 -2723.2833 874.9594 11950.115
1000 310.3971 -3181.9675 -2720.2753 842.81184 11950.115
Loop time of 0.279202 on 1 procs for 1000 steps with 500 atoms
Performance: 309.453 ns/day, 0.078 hours/ns, 3581.633 timesteps/s
99.1% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.24196 | 0.24196 | 0.24196 | 0.0 | 86.66
Bond | 6.628e-05 | 6.628e-05 | 6.628e-05 | 0.0 | 0.02
Neigh | 0.0043204 | 0.0043204 | 0.0043204 | 0.0 | 1.55
Comm | 0.010242 | 0.010242 | 0.010242 | 0.0 | 3.67
Output | 0.0012252 | 0.0012252 | 0.0012252 | 0.0 | 0.44
Modify | 0.017572 | 0.017572 | 0.017572 | 0.0 | 6.29
Other | | 0.003811 | | | 1.37
Nlocal: 500 ave 500 max 500 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1610 ave 1610 max 1610 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 14767 ave 14767 max 14767 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 14767
Ave neighs/atom = 29.534
Ave special neighs/atom = 0
Neighbor list builds = 5
Dangerous builds = 0
#write_data lj-out.data
Total wall time: 0:00:00

View File

@ -0,0 +1,173 @@
LAMMPS (9 Nov 2016)
using 1 OpenMP thread(s) per MPI task
# LJ particles
variable T0 equal 300.0
variable press equal 0.0
units real
atom_style full
pair_style lj/cut 5.0
read_data "lj.data"
orthogonal box = (1.06874 1.06874 1.06874) to (23.9313 23.9313 23.9313)
1 by 2 by 2 MPI processor grid
reading atoms ...
500 atoms
reading velocities ...
500 velocities
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
thermo 10
thermo_style custom step temp pe etotal press vol
timestep 1.0
fix fxnvt all nvt temp ${T0} ${T0} 1000.0
fix fxnvt all nvt temp 300 ${T0} 1000.0
fix fxnvt all nvt temp 300 300 1000.0
fix fxgREM all grem 400 -.01 -30000 fxnvt
run 1000
Neighbor list info ...
1 neighbor list requests
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 7
ghost atom cutoff = 7
binsize = 3.5 -> bins = 7 7 7
Memory usage per processor = 5.34276 Mbytes
Step Temp PotEng TotEng Press Volume
0 305.69499 -3177.6423 -2722.9442 883.58369 11950.115
10 312.30121 -3182.2257 -2717.7013 793.47811 11950.115
20 314.94553 -3186.4559 -2717.9983 738.74091 11950.115
30 312.22861 -3183.7638 -2719.3474 797.47978 11950.115
40 305.93987 -3180.7079 -2725.6455 881.30806 11950.115
50 300.42132 -3176.5828 -2729.7288 967.92042 11950.115
60 299.16487 -3174.1921 -2729.2071 1004.247 11950.115
70 301.65565 -3176.0891 -2727.3992 962.58134 11950.115
80 304.77334 -3178.2663 -2724.939 907.8946 11950.115
90 305.59929 -3178.9472 -2724.3914 879.91629 11950.115
100 303.79263 -3177.5103 -2725.6418 892.67631 11950.115
110 300.85863 -3175.4711 -2727.9667 923.44924 11950.115
120 298.69083 -3175.6415 -2731.3615 931.87518 11950.115
130 298.38415 -3176.3706 -2732.5468 928.88286 11950.115
140 299.99129 -3177.6935 -2731.4792 914.36783 11950.115
150 301.83869 -3178.121 -2729.1588 915.01407 11950.115
160 301.18834 -3177.3117 -2729.3169 947.45228 11950.115
170 296.99406 -3172.7363 -2730.9801 1042.6928 11950.115
180 291.25952 -3168.3407 -2735.1142 1144.5436 11950.115
190 287.1178 -3164.9847 -2737.9187 1223.4003 11950.115
200 286.552 -3165.2799 -2739.0555 1235.6703 11950.115
210 289.81033 -3167.8062 -2736.7353 1194.6672 11950.115
220 296.09616 -3171.7847 -2731.3641 1115.8799 11950.115
230 303.79176 -3176.2893 -2724.4221 1024.6471 11950.115
240 309.9273 -3180.9591 -2719.9657 945.55045 11950.115
250 312.0159 -3182.3307 -2718.2306 934.36956 11950.115
260 309.63264 -3181.0452 -2720.4901 986.77385 11950.115
270 304.38172 -3176.568 -2723.8233 1097.264 11950.115
280 298.64188 -3174.1384 -2729.9313 1186.2239 11950.115
290 294.76686 -3170.9562 -2732.5128 1264.247 11950.115
300 294.00805 -3169.8091 -2732.4944 1287.4001 11950.115
310 296.41801 -3172.834 -2731.9347 1229.5624 11950.115
320 301.38477 -3175.9514 -2727.6644 1140.8664 11950.115
330 307.86584 -3181.2442 -2723.3171 1007.1545 11950.115
340 314.7103 -3185.9891 -2717.8814 871.74528 11950.115
350 320.53954 -3187.8385 -2711.0602 776.85994 11950.115
360 323.49505 -3188.9927 -2707.8184 716.58062 11950.115
370 321.62077 -3187.1246 -2708.7381 731.01909 11950.115
380 314.39049 -3181.4931 -2713.8611 831.21057 11950.115
390 303.57079 -3174.1804 -2722.6419 978.62645 11950.115
400 293.42165 -3167.0452 -2730.6027 1122.3558 11950.115
410 288.46838 -3165.4071 -2736.3322 1171.8087 11950.115
420 290.30766 -3168.5988 -2736.7882 1122.5413 11950.115
430 296.34338 -3173.0824 -2732.2941 1030.2769 11950.115
440 301.92394 -3175.4307 -2726.3417 964.25387 11950.115
450 303.76745 -3176.9122 -2725.0811 934.49176 11950.115
460 301.72985 -3174.2821 -2725.4818 979.07605 11950.115
470 298.93736 -3173.0014 -2728.3548 1020.0482 11950.115
480 298.80912 -3172.803 -2728.3471 1036.6531 11950.115
490 302.72217 -3175.3764 -2725.1001 997.71146 11950.115
500 309.11393 -3179.3088 -2719.5253 925.81108 11950.115
510 314.37612 -3183.2961 -2715.6855 856.23748 11950.115
520 315.72767 -3183.547 -2713.926 847.70543 11950.115
530 313.34173 -3182.2695 -2716.1974 877.30842 11950.115
540 309.84312 -3178.9553 -2718.0871 936.69244 11950.115
550 308.3251 -3177.8582 -2719.248 963.93032 11950.115
560 308.79192 -3176.4834 -2717.1788 989.67643 11950.115
570 307.57194 -3175.8464 -2718.3565 1021.0494 11950.115
580 301.8035 -3173.1582 -2724.2483 1102.4893 11950.115
590 292.43425 -3165.751 -2730.7772 1254.7815 11950.115
600 283.62905 -3159.8987 -2738.022 1381.0608 11950.115
610 279.90122 -3157.49 -2741.1581 1431.0028 11950.115
620 283.4582 -3160.756 -2739.1334 1367.7385 11950.115
630 292.58866 -3166.3469 -2731.1435 1241.1194 11950.115
640 302.86585 -3173.4778 -2722.9878 1089.7342 11950.115
650 309.89252 -3179.4078 -2718.4662 972.6359 11950.115
660 312.19165 -3182.7754 -2718.414 916.62037 11950.115
670 311.2287 -3181.9102 -2718.9811 933.79804 11950.115
680 308.57852 -3180.7312 -2721.7441 969.24936 11950.115
690 304.71609 -3176.9196 -2723.6775 1040.2699 11950.115
700 300.31995 -3175.5245 -2728.8213 1082.845 11950.115
710 296.43537 -3173.0166 -2732.0915 1127.4487 11950.115
720 293.86692 -3172.1582 -2735.0535 1135.0215 11950.115
730 293.35611 -3170.9335 -2734.5885 1122.9143 11950.115
740 295.44861 -3172.9862 -2733.5288 1050.995 11950.115
750 299.82732 -3175.3467 -2729.3763 958.31462 11950.115
760 305.37987 -3178.216 -2723.9866 854.1946 11950.115
770 310.75394 -3182.8127 -2720.5898 737.72668 11950.115
780 314.81395 -3183.7905 -2715.5286 679.74198 11950.115
790 316.63339 -3185.8028 -2714.8346 638.48871 11950.115
800 315.2894 -3186.9345 -2717.9654 641.53256 11950.115
810 310.54289 -3183.7383 -2721.8293 728.51241 11950.115
820 303.31439 -3178.7897 -2727.6326 864.45674 11950.115
830 295.46125 -3175.5387 -2736.0625 997.72969 11950.115
840 288.802 -3169.6502 -2740.0791 1160.6622 11950.115
850 283.94785 -3165.8605 -2743.5096 1289.55 11950.115
860 280.17501 -3163.0381 -2746.299 1392.8854 11950.115
870 277.2456 -3159.8429 -2747.4611 1481.3899 11950.115
880 275.93123 -3158.3584 -2747.9316 1523.5374 11950.115
890 277.0215 -3159.2285 -2747.18 1506.1558 11950.115
900 280.96237 -3162.483 -2744.5728 1428.4183 11950.115
910 287.37962 -3167.6183 -2740.1628 1303.0268 11950.115
920 294.56731 -3171.6765 -2733.5299 1177.748 11950.115
930 300.63273 -3174.0842 -2726.9158 1078.7393 11950.115
940 304.10943 -3175.9847 -2723.645 1007.7154 11950.115
950 304.64845 -3176.6263 -2723.4848 976.37917 11950.115
960 303.36343 -3176.4694 -2725.2393 971.40749 11950.115
970 302.57138 -3177.5541 -2727.5021 954.01115 11950.115
980 304.2593 -3179.2101 -2726.6475 919.74949 11950.115
990 307.69959 -3180.9631 -2723.2833 874.9594 11950.115
1000 310.3971 -3181.9675 -2720.2753 842.81184 11950.115
Loop time of 0.133894 on 4 procs for 1000 steps with 500 atoms
Performance: 645.285 ns/day, 0.037 hours/ns, 7468.580 timesteps/s
98.8% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.065271 | 0.071043 | 0.07818 | 1.9 | 53.06
Bond | 5.6505e-05 | 6.5565e-05 | 7.7724e-05 | 0.1 | 0.05
Neigh | 0.0011396 | 0.0012607 | 0.0013669 | 0.2 | 0.94
Comm | 0.033866 | 0.040269 | 0.045386 | 2.6 | 30.08
Output | 0.0019252 | 0.0020776 | 0.0023642 | 0.4 | 1.55
Modify | 0.012141 | 0.013629 | 0.01486 | 0.9 | 10.18
Other | | 0.005549 | | | 4.14
Nlocal: 125 ave 127 max 123 min
Histogram: 1 0 1 0 0 0 0 1 0 1
Nghost: 871.25 ave 882 max 863 min
Histogram: 2 0 0 0 0 0 1 0 0 1
Neighs: 3691.75 ave 3808 max 3563 min
Histogram: 1 0 0 0 1 0 1 0 0 1
Total # of neighbors = 14767
Ave neighs/atom = 29.534
Ave special neighs/atom = 0
Neighbor list builds = 5
Dangerous builds = 0
#write_data lj-out.data
Total wall time: 0:00:00

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

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,31 @@
# LJ particles
variable lambda world 900 910 920 930
variable rep world 0 1 2 3
#variable walker world 0 1 3 2
variable T0 equal 300.0
variable press equal 0.0
units real
atom_style full
pair_style lj/cut 5.0
# LJ particles
log ${rep}/log.lammps.${rep}
print "This is replica: ${rep}"
read_data ${rep}/lj.data
#dump dump all xyz 1000 ${rep}/dump.xyz
thermo 10
thermo_style custom step temp pe etotal press vol
timestep 1.0
fix fxnpt all npt temp ${T0} ${T0} 1000.0 iso ${press} ${press} 10000.0
fix fxgREM all grem ${lambda} -.03 -30000 fxnpt
thermo_modify press fxgREM_press
temper/grem 10000 100 ${lambda} fxgREM fxnpt 10294 98392 #${walker}
#write_data ${rep}/lj-out.data

6
src/.gitignore vendored
View File

@ -205,6 +205,8 @@
/compute_pe_tally.h
/compute_plasticity_atom.cpp
/compute_plasticity_atom.h
/compute_pressure_grem.cpp
/compute_pressure_grem.h
/compute_rigid_local.cpp
/compute_rigid_local.h
/compute_spec_atom.cpp
@ -343,6 +345,8 @@
/fix_gle.h
/fix_gpu.cpp
/fix_gpu.h
/fix_grem.cpp
/fix_grem.h
/fix_imd.cpp
/fix_imd.h
/fix_ipi.cpp
@ -873,6 +877,8 @@
/tad.h
/temper.cpp
/temper.h
/temper_grem.cpp
/temper_grem.h
/thr_data.cpp
/thr_data.h
/verlet_split.cpp

View File

@ -0,0 +1,156 @@
/* ----------------------------------------------------------------------
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 <mpi.h>
#include <string.h>
#include <stdlib.h>
#include "compute_pressure_grem.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "modify.h"
#include "fix.h"
#include "force.h"
#include "pair.h"
#include "kspace.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ----------------------------------------------------------------------
Last argument is the id of the gREM fix
------------------------------------------------------------------------- */
ComputePressureGrem::ComputePressureGrem(LAMMPS *lmp, int narg, char **arg) :
ComputePressure(lmp, narg-1, arg)
{
int len = strlen(arg[narg-1])+1;
fix_grem = new char[len];
strcpy(fix_grem,arg[narg-1]);
}
/* ---------------------------------------------------------------------- */
ComputePressureGrem::~ComputePressureGrem()
{
delete [] fix_grem;
}
/* ---------------------------------------------------------------------- */
void ComputePressureGrem::init()
{
ComputePressure::init();
// Initialize hook to gREM fix
int ifix = modify->find_fix(fix_grem);
if (ifix < 0)
error->all(FLERR,"Fix grem ID for compute pressure/grem does not exist");
int dim;
scale_grem = (double *)modify->fix[ifix]->extract("scale_grem",dim);
if (scale_grem == NULL || dim != 0)
error->all(FLERR,"Cannot extract gREM scale factor from fix grem");
}
/* ----------------------------------------------------------------------
compute total pressure, averaged over Pxx, Pyy, Pzz
------------------------------------------------------------------------- */
double ComputePressureGrem::compute_scalar()
{
invoked_scalar = update->ntimestep;
if (update->vflag_global != invoked_scalar)
error->all(FLERR,"Virial was not tallied on needed timestep");
// invoke temperature if it hasn't been already
double t;
if (keflag) {
if (temperature->invoked_scalar != update->ntimestep)
t = temperature->compute_scalar() / (*scale_grem);
else t = temperature->scalar / (*scale_grem);
}
if (dimension == 3) {
inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
virial_compute(3,3);
if (keflag)
scalar = (temperature->dof * boltz * t +
virial[0] + virial[1] + virial[2]) / 3.0 * inv_volume * nktv2p;
else
scalar = (virial[0] + virial[1] + virial[2]) / 3.0 * inv_volume * nktv2p;
} else {
inv_volume = 1.0 / (domain->xprd * domain->yprd);
virial_compute(2,2);
if (keflag)
scalar = (temperature->dof * boltz * t +
virial[0] + virial[1]) / 2.0 * inv_volume * nktv2p;
else
scalar = (virial[0] + virial[1]) / 2.0 * inv_volume * nktv2p;
}
return scalar;
}
/* ----------------------------------------------------------------------
compute pressure tensor
assume KE tensor has already been computed
------------------------------------------------------------------------- */
void ComputePressureGrem::compute_vector()
{
invoked_vector = update->ntimestep;
if (update->vflag_global != invoked_vector)
error->all(FLERR,"Virial was not tallied on needed timestep");
if (force->kspace && kspace_virial && force->kspace->scalar_pressure_flag)
error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' for "
"tensor components with kspace_style msm");
// invoke temperature if it hasn't been already
double ke_tensor[6];
if (keflag) {
if (temperature->invoked_vector != update->ntimestep)
temperature->compute_vector();
for (int i = 0; i < 6; ++i)
ke_tensor[i] = temperature->vector[i] / (*scale_grem);
}
if (dimension == 3) {
inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
virial_compute(6,3);
if (keflag) {
for (int i = 0; i < 6; i++)
vector[i] = (ke_tensor[i] + virial[i]) * inv_volume * nktv2p;
} else
for (int i = 0; i < 6; i++)
vector[i] = virial[i] * inv_volume * nktv2p;
} else {
inv_volume = 1.0 / (domain->xprd * domain->yprd);
virial_compute(4,2);
if (keflag) {
vector[0] = (ke_tensor[0] + virial[0]) * inv_volume * nktv2p;
vector[1] = (ke_tensor[1] + virial[1]) * inv_volume * nktv2p;
vector[3] = (ke_tensor[3] + virial[3]) * inv_volume * nktv2p;
vector[2] = vector[4] = vector[5] = 0.0;
} else {
vector[0] = virial[0] * inv_volume * nktv2p;
vector[1] = virial[1] * inv_volume * nktv2p;
vector[3] = virial[3] * inv_volume * nktv2p;
vector[2] = vector[4] = vector[5] = 0.0;
}
}
}

View File

@ -0,0 +1,91 @@
/* -*- 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.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(pressure/grem,ComputePressureGrem)
#else
#ifndef LMP_COMPUTE_PRESSURE_GREM_H
#define LMP_COMPUTE_PRESSURE_GREM_H
#include "compute_pressure.h"
namespace LAMMPS_NS {
class ComputePressureGrem : public ComputePressure {
public:
ComputePressureGrem(class LAMMPS *, int, char **);
virtual ~ComputePressureGrem();
virtual void init();
virtual double compute_scalar();
virtual void compute_vector();
protected:
// Access to gREM fix scale factor
char *fix_grem;
double *scale_grem;
};
}
#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.
E: Compute pressure must use group all
Virial contributions computed by potentials (pair, bond, etc) are
computed on all atoms.
E: Could not find compute pressure temperature ID
The compute ID for calculating temperature does not exist.
E: Compute pressure temperature ID does not compute temperature
The compute ID assigned to a pressure computation must compute
temperature.
E: Compute pressure requires temperature ID to include kinetic energy
The keflag cannot be used unless a temperature compute is provided.
E: Virial was not tallied on needed timestep
You are using a thermo keyword that requires potentials to
have tallied the virial, but they didn't on this timestep. See the
variable doc page for ideas on how to make this work.
E: Must use 'kspace_modify pressure/scalar no' for tensor components with kspace_style msm
Otherwise MSM will compute only a scalar pressure. See the kspace_modify
command for details on this setting.
E: Fix grem ID for compute pressure/grem does not exist
Compute pressure/grem was passed an invalid fix id
E: Cannot extract gREM scale factor from fix grem
The fix id passed to compute pressure/grem refers to an incompatible fix
*/

296
src/USER-MISC/fix_grem.cpp Normal file
View File

@ -0,0 +1,296 @@
/* ----------------------------------------------------------------------
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.
Force scaling fix for gREM.
Cite: http://dx.doi.org/10.1063/1.3432176
Cite: http://dx.doi.org/10.1021/acs.jpcb.5b07614
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Edyta Malolepsza (Broad Institute)
Contributing author: David Stelter (Boston University)
Contributing author: Tom Keyes (Boston University)
------------------------------------------------------------------------- */
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include "comm.h"
#include "fix_grem.h"
#include "atom.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "domain.h"
#include "input.h"
#include "compute.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum{NONE,CONSTANT,EQUAL,ATOM};
/* ---------------------------------------------------------------------- */
FixGrem::FixGrem(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 7) error->all(FLERR,"Illegal fix grem command");
scalar_flag = 1;
vector_flag = 1;
size_vector = 3;
global_freq = 1;
extscalar = 1;
extvector = 1;
scale_grem = 1.0;
// tbath - temp of bath, the same as defined in thermostat
lambda = force->numeric(FLERR,arg[3]);
eta = force->numeric(FLERR,arg[4]);
h0 = force->numeric(FLERR,arg[5]);
int n = strlen(arg[6])+1;
id_nh = new char[n];
strcpy(id_nh,arg[6]);
// create a new compute temp style
// id = fix-ID + temp
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
newarg = new char*[5];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure/grem";
newarg[3] = id_temp;
newarg[4] = id;
modify->add_compute(5,newarg);
delete [] newarg;
// create a new compute ke style
// id = fix-ID + ke
n = strlen(id) + 8;
id_ke = new char[n];
strcpy(id_ke,id);
strcat(id_ke,"_ke");
newarg = new char*[3];
newarg[0] = id_ke;
newarg[1] = (char *) "all";
newarg[2] = (char *) "ke";
modify->add_compute(3,newarg);
delete [] newarg;
// create a new compute pe style
// id = fix-ID + pe
n = strlen(id) + 9;
id_pe = new char[n];
strcpy(id_pe,id);
strcat(id_pe,"_pe");
newarg = new char*[3];
newarg[0] = id_pe;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pe";
modify->add_compute(3,newarg);
delete [] newarg;
int ifix = modify->find_fix(id_nh);
if (ifix < 0)
error->all(FLERR,"Fix id for nvt or npt fix does not exist");
Fix *nh = modify->fix[ifix];
pressflag = 0;
int *p_flag = (int *)nh->extract("p_flag",ifix);
if ((p_flag == NULL) || (ifix != 1) || (p_flag[0] == 0)
|| (p_flag[1] == 0) || (p_flag[2] == 0)) {
pressflag = 0;
} else if ((p_flag[0] == 1) && (p_flag[1] == 1)
&& (p_flag[2] == 1) && (ifix == 1)) {
pressflag = 1;
char *modargs[2];
modargs[0] = (char *) "press";
modargs[1] = id_press;
nh->modify_param(2,modargs);
}
}
/* ---------------------------------------------------------------------- */
FixGrem::~FixGrem()
{
// delete temperature, pressure and energies if fix created them
modify->delete_compute(id_temp);
modify->delete_compute(id_press);
modify->delete_compute(id_ke);
modify->delete_compute(id_pe);
delete [] id_temp;
delete [] id_press;
delete [] id_ke;
delete [] id_pe;
delete [] id_nh;
}
/* ---------------------------------------------------------------------- */
int FixGrem::setmask()
{
int mask = 0;
mask |= POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixGrem::init()
{
if (domain->triclinic)
error->all(FLERR,"Triclinic cells are not supported");
// set temperature and pressure ptrs
int icompute = modify->find_compute(id_temp);
if (icompute < 0)
error->all(FLERR,"Temperature compute ID for fix grem does not exist");
temperature = modify->compute[icompute];
icompute = modify->find_compute(id_ke);
if (icompute < 0)
error->all(FLERR,"KE compute ID for fix grem does not exist");
ke = modify->compute[icompute];
icompute = modify->find_compute(id_pe);
if (icompute < 0)
error->all(FLERR,"PE compute ID for fix grem does not exist");
pe = modify->compute[icompute];
int ifix = modify->find_fix(id_nh);
if (ifix < 0)
error->all(FLERR,"Fix id for nvt or npt fix does not exist");
Fix *nh = modify->fix[ifix];
double *t_start = (double *)nh->extract("t_start",ifix);
double *t_stop = (double *)nh->extract("t_stop",ifix);
if ((t_start != NULL) && (t_stop != NULL) && (ifix == 0)) {
tbath = *t_start;
if (*t_start != *t_stop)
error->all(FLERR,"Thermostat temperature ramp not allowed");
} else
error->all(FLERR,"Problem extracting target temperature from fix nvt or npt");
pressref = 0.0;
if (pressflag) {
int *p_flag = (int *)nh->extract("p_flag",ifix);
double *p_start = (double *) nh->extract("p_start",ifix);
double *p_stop = (double *) nh->extract("p_stop",ifix);
if ((p_flag != NULL) && (p_start != NULL) && (p_stop != NULL)
&& (ifix == 1)) {
ifix = 0;
pressref = p_start[0];
if ((p_start[0] != p_stop[0]) || (p_flag[0] != 1)) ++ ifix;
if ((p_start[1] != p_stop[1]) || (p_flag[0] != 1)) ++ ifix;
if ((p_start[2] != p_stop[2]) || (p_flag[0] != 1)) ++ ifix;
if ((p_start[0] != p_start[1]) || (p_start[1] != p_start[2])) ++ifix;
if ((p_flag[3] != 0) || (p_flag[4] != 0) || (p_flag[5] != 0)) ++ifix;
if (ifix > 0)
error->all(FLERR,"Unsupported pressure settings in fix npt");
} else
error->all(FLERR,"Problem extracting target pressure from fix npt");
}
}
/* ---------------------------------------------------------------------- */
void FixGrem::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
if (strstr(update->integrate_style,"respa"))
error->all(FLERR,"Run style 'respa' is not supported");
}
/* ---------------------------------------------------------------------- */
void FixGrem::min_setup(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixGrem::post_force(int vflag)
{
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double tmpvolume = domain->xprd * domain->yprd * domain->zprd;
double tmppe = pe->compute_scalar();
// potential energy
double tmpenthalpy = tmppe+pressref*tmpvolume/(force->nktv2p);
double teffective = lambda+eta*(tmpenthalpy-h0);
scale_grem = tbath/teffective;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
f[i][0] *= scale_grem;
f[i][1] *= scale_grem;
f[i][2] *= scale_grem;
}
pe->addstep(update->ntimestep+1);
}
/* ----------------------------------------------------------------------
extract scale factor
------------------------------------------------------------------------- */
void *FixGrem::extract(const char *str, int &dim)
{
dim=0;
if (strcmp(str,"scale_grem") == 0) {
return &scale_grem;
}
return NULL;
}

83
src/USER-MISC/fix_grem.h Normal file
View File

@ -0,0 +1,83 @@
/* -*- 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.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(grem,FixGrem)
#else
#ifndef LMP_FIX_GREM_H
#define LMP_FIX_GREM_H
#include "fix.h"
namespace LAMMPS_NS {
class FixGrem : public Fix {
public:
FixGrem(class LAMMPS *, int, char **);
~FixGrem();
int setmask();
void init();
void setup(int);
void min_setup(int);
void post_force(int);
void *extract(const char *, int &);
double scale_grem,lambda,eta,h0;
int pressflag;
private:
double tbath,pressref;
protected:
char *id_temp,*id_press,*id_ke,*id_pe,*id_nh;
class Compute *temperature,*pressure,*ke,*pe;
};
}
#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.
E: Region ID for fix grem does not exist
Self-explanatory.
E: Variable name for fix grem does not exist
Self-explanatory.
E: Variable for fix grem is invalid style
Self-explanatory.
E: Cannot use variable energy with constant force in fix grem
This is because for constant force, LAMMPS can compute the change
in energy directly.
E: Must use variable energy with fix grem
Must define an energy vartiable when applyting a dynamic
force during minimization.
*/

View File

@ -0,0 +1,378 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Mark Sears (SNL)
Modified for gREM by David Stelter (BU)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "temper_grem.h"
#include "fix_grem.h"
#include "universe.h"
#include "domain.h"
#include "atom.h"
#include "update.h"
#include "integrate.h"
#include "modify.h"
#include "compute.h"
#include "force.h"
#include "output.h"
#include "thermo.h"
#include "fix.h"
#include "random_park.h"
#include "finish.h"
#include "timer.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
//#define TEMPER_DEBUG 1
/* ---------------------------------------------------------------------- */
TemperGrem::TemperGrem(LAMMPS *lmp) : Pointers(lmp) {}
/* ---------------------------------------------------------------------- */
TemperGrem::~TemperGrem()
{
MPI_Comm_free(&roots);
if (ranswap) delete ranswap;
delete ranboltz;
delete [] set_lambda;
delete [] lambda2world;
delete [] world2lambda;
delete [] world2root;
delete [] id_nh;
}
/* ----------------------------------------------------------------------
perform tempering with inter-world swaps
------------------------------------------------------------------------- */
void TemperGrem::command(int narg, char **arg)
{
if (universe->nworlds == 1)
error->all(FLERR,"Must have more than one processor partition to temper");
if (domain->box_exist == 0)
error->all(FLERR,"Temper/gREM command before simulation box is defined");
if (narg != 7 && narg != 8)
error->universe_all(FLERR,"Illegal temper command");
int nsteps = force->inumeric(FLERR,arg[0]);
nevery = force->inumeric(FLERR,arg[1]);
double lambda = force->numeric(FLERR,arg[2]);
// Get and check if gREM fix exists
for (whichfix = 0; whichfix < modify->nfix; whichfix++)
if (strcmp(arg[3],modify->fix[whichfix]->id) == 0) break;
if (whichfix == modify->nfix)
error->universe_all(FLERR,"Tempering fix ID is not defined");
fix_grem = (FixGrem*)(modify->fix[whichfix]);
// Check input values lambdas should be equal, assign other gREM values
if (lambda != fix_grem->lambda)
error->universe_all(FLERR,"Lambda from tempering and fix in the same world"
" must be the same");
double eta = fix_grem->eta;
double h0 = fix_grem->h0;
double pressref = 0;
// Get and check for nh fix
int n = strlen(arg[4])+1;
id_nh = new char[n];
strcpy(id_nh,arg[4]);
int ifix = modify->find_fix(id_nh);
if (ifix < 0)
error->all(FLERR,"Fix id for nvt or npt fix does not exist");
Fix *nh = modify->fix[ifix];
// get result from nvt vs npt check from fix_grem
int pressflag = fix_grem->pressflag;
// fix_grem does all the checking...
if (pressflag) {
double *p_start = (double *) nh->extract("p_start",ifix);
pressref = p_start[0];
}
seed_swap = force->inumeric(FLERR,arg[5]);
seed_boltz = force->inumeric(FLERR,arg[6]);
my_set_lambda = universe->iworld;
if (narg == 8) my_set_lambda = force->inumeric(FLERR,arg[7]);
// swap frequency must evenly divide total # of timesteps
if (nevery <= 0)
error->universe_all(FLERR,"Invalid frequency in temper command");
nswaps = nsteps/nevery;
if (nswaps*nevery != nsteps)
error->universe_all(FLERR,"Non integer # of swaps in temper command");
// Must be used with fix_grem
if (strcmp(modify->fix[whichfix]->style,"grem") != 0)
error->universe_all(FLERR,"Tempering temperature fix is not supported");
// setup for long tempering run
update->whichflag = 1;
update->nsteps = nsteps;
update->beginstep = update->firststep = update->ntimestep;
update->endstep = update->laststep = update->firststep + nsteps;
if (update->laststep < 0)
error->all(FLERR,"Too many timesteps");
lmp->init();
// local storage
me_universe = universe->me;
MPI_Comm_rank(world,&me);
nworlds = universe->nworlds;
iworld = universe->iworld;
boltz = force->boltz;
// pe_compute = ptr to thermo_pe compute
// notify compute it will be called at first swap
int id = modify->find_compute("thermo_pe");
if (id < 0) error->all(FLERR,"Tempering could not find thermo_pe compute");
Compute *pe_compute = modify->compute[id];
pe_compute->addstep(update->ntimestep + nevery);
// create MPI communicator for root proc from each world
int color;
if (me == 0) color = 0;
else color = 1;
MPI_Comm_split(universe->uworld,color,0,&roots);
// RNGs for swaps and Boltzmann test
// warm up Boltzmann RNG
if (seed_swap) ranswap = new RanPark(lmp,seed_swap);
else ranswap = NULL;
ranboltz = new RanPark(lmp,seed_boltz + me_universe);
for (int i = 0; i < 100; i++) ranboltz->uniform();
// world2root[i] = global proc that is root proc of world i
world2root = new int[nworlds];
if (me == 0)
MPI_Allgather(&me_universe,1,MPI_INT,world2root,1,MPI_INT,roots);
MPI_Bcast(world2root,nworlds,MPI_INT,0,world);
// create static list of set lambdas
// allgather tempering arg "lambda" across root procs
// bcast from each root to other procs in world
set_lambda = new double[nworlds];
if (me == 0) MPI_Allgather(&lambda,1,MPI_DOUBLE,set_lambda,1,MPI_DOUBLE,roots);
MPI_Bcast(set_lambda,nworlds,MPI_DOUBLE,0,world);
// create world2lambda only on root procs from my_set_lambda
// create lambda2world on root procs from world2lambda,
// then bcast to all procs within world
world2lambda = new int[nworlds];
lambda2world = new int[nworlds];
if (me == 0) {
MPI_Allgather(&my_set_lambda,1,MPI_INT,world2lambda,1,MPI_INT,roots);
for (int i = 0; i < nworlds; i++) lambda2world[world2lambda[i]] = i;
}
MPI_Bcast(lambda2world,nworlds,MPI_INT,0,world);
// if restarting tempering, reset lambda target of Fix to current my_set_lambda
if (narg == 8) {
double new_lambda = set_lambda[my_set_lambda];
fix_grem->lambda = new_lambda;
}
// setup tempering runs
int i,which,partner,swap,partner_set_lambda,partner_world;
double pe,weight,weight_partner,weight_cross, weight_cross_partner;
double volume,enth,new_lambda,boltz_factor;
if (me_universe == 0 && universe->uscreen)
fprintf(universe->uscreen,"Setting up tempering ...\n");
update->integrate->setup();
if (me_universe == 0) {
if (universe->uscreen) {
fprintf(universe->uscreen,"Step");
for (int i = 0; i < nworlds; i++)
fprintf(universe->uscreen," T%d",i);
fprintf(universe->uscreen,"\n");
}
if (universe->ulogfile) {
fprintf(universe->ulogfile,"Step");
for (int i = 0; i < nworlds; i++)
fprintf(universe->ulogfile," T%d",i);
fprintf(universe->ulogfile,"\n");
}
print_status();
}
timer->init();
timer->barrier_start();
for (int iswap = 0; iswap < nswaps; iswap++) {
// run for nevery timesteps
update->integrate->run(nevery);
// compute PE
// notify compute it will be called at next swap
pe = pe_compute->compute_scalar();
pe_compute->addstep(update->ntimestep + nevery);
// which = which of 2 kinds of swaps to do (0,1)
if (!ranswap) which = iswap % 2;
else if (ranswap->uniform() < 0.5) which = 0;
else which = 1;
// partner_set_lambda = which set lambda I am partnering with for this swap
if (which == 0) {
if (my_set_lambda % 2 == 0) partner_set_lambda = my_set_lambda + 1;
else partner_set_lambda = my_set_lambda - 1;
} else {
if (my_set_lambda % 2 == 1) partner_set_lambda = my_set_lambda + 1;
else partner_set_lambda = my_set_lambda - 1;
}
// partner = proc ID to swap with
// if partner = -1, then I am not a proc that swaps
partner = -1;
if (me == 0 && partner_set_lambda >= 0 && partner_set_lambda < nworlds) {
partner_world = lambda2world[partner_set_lambda];
partner = world2root[partner_world];
}
// compute weights
volume = domain->xprd * domain->yprd * domain->zprd;
enth = pe + (pressref * volume);
weight = log(set_lambda[my_set_lambda] + (eta*(enth + h0)));
weight_cross = log(set_lambda[partner_set_lambda] + (eta*(enth + h0)));
// swap with a partner, only root procs in each world participate
// hi proc sends PE to low proc
// lo proc make Boltzmann decision on whether to swap
// lo proc communicates decision back to hi proc
swap = 0;
if (partner != -1) {
if (me_universe > partner) {
MPI_Send(&weight,1,MPI_DOUBLE,partner,0,universe->uworld);
MPI_Send(&weight_cross,1,MPI_DOUBLE,partner,0,universe->uworld);
}
else {
MPI_Recv(&weight_partner,1,MPI_DOUBLE,partner,0,universe->uworld,MPI_STATUS_IGNORE);
MPI_Recv(&weight_cross_partner,1,MPI_DOUBLE,partner,0,universe->uworld,MPI_STATUS_IGNORE);
}
if (me_universe < partner) {
boltz_factor = (weight - weight_partner + weight_cross - weight_cross_partner) *
(1 / (boltz * eta));
if (boltz_factor >= 0.0) swap = 1;
else if (ranboltz->uniform() < exp(boltz_factor)) swap = 1;
}
if (me_universe < partner)
MPI_Send(&swap,1,MPI_INT,partner,0,universe->uworld);
else
MPI_Recv(&swap,1,MPI_INT,partner,0,universe->uworld,MPI_STATUS_IGNORE);
#ifdef TEMPER_DEBUG
if (me_universe < partner)
printf("SWAP %d & %d: yes = %d,Ts = %d %d, PEs = %g %g, Bz = %g %g\n",
me_universe,partner,swap,my_set_lambda,partner_set_lambda,
weight,weight_partner,boltz_factor,exp(boltz_factor));
#endif
}
// bcast swap result to other procs in my world
MPI_Bcast(&swap,1,MPI_INT,0,world);
// if my world swapped, all procs in world reset temp target of Fix
if (swap) {
new_lambda = set_lambda[partner_set_lambda];
fix_grem->lambda = new_lambda;
}
// update my_set_lambda and lambda2world on every proc
// root procs update their value if swap took place
// allgather across root procs
// bcast within my world
if (swap) my_set_lambda = partner_set_lambda;
if (me == 0) {
MPI_Allgather(&my_set_lambda,1,MPI_INT,world2lambda,1,MPI_INT,roots);
for (i = 0; i < nworlds; i++) lambda2world[world2lambda[i]] = i;
}
MPI_Bcast(lambda2world,nworlds,MPI_INT,0,world);
// print out current swap status
if (me_universe == 0) print_status();
}
timer->barrier_stop();
update->integrate->cleanup();
Finish finish(lmp);
finish.end(1);
update->whichflag = 0;
update->firststep = update->laststep = 0;
update->beginstep = update->endstep = 0;
}
/* ----------------------------------------------------------------------
proc 0 prints current tempering status
------------------------------------------------------------------------- */
void TemperGrem::print_status()
{
if (universe->uscreen) {
fprintf(universe->uscreen,BIGINT_FORMAT,update->ntimestep);
for (int i = 0; i < nworlds; i++)
fprintf(universe->uscreen," %d",world2lambda[i]);
fprintf(universe->uscreen,"\n");
}
if (universe->ulogfile) {
fprintf(universe->ulogfile,BIGINT_FORMAT,update->ntimestep);
for (int i = 0; i < nworlds; i++)
fprintf(universe->ulogfile," %d",world2lambda[i]);
fprintf(universe->ulogfile,"\n");
fflush(universe->ulogfile);
}
}

111
src/USER-MISC/temper_grem.h Normal file
View File

@ -0,0 +1,111 @@
/* -*- 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.
------------------------------------------------------------------------- */
#ifdef COMMAND_CLASS
CommandStyle(temper/grem,TemperGrem)
#else
#ifndef LMP_TEMPER_GREM_H
#define LMP_TEMPER_GREM_H
#include "pointers.h"
namespace LAMMPS_NS {
class TemperGrem : protected Pointers {
public:
TemperGrem(class LAMMPS *);
~TemperGrem();
void command(int, char **);
private:
int me,me_universe; // my proc ID in world and universe
int iworld,nworlds; // world info
double boltz; // copy from output->boltz
MPI_Comm roots; // MPI comm with 1 root proc from each world
class RanPark *ranswap,*ranboltz; // RNGs for swapping and Boltz factor
int nevery; // # of timesteps between swaps
int nswaps; // # of tempering swaps to perform
int seed_swap; // 0 = toggle swaps, n = RNG for swap direction
int seed_boltz; // seed for Boltz factor comparison
int whichfix; // index of temperature fix to use
int fixstyle; // what kind of temperature fix is used
int my_set_lambda; // which set lambda I am simulating
double *set_lambda; // static list of replica set lambdas
int *lambda2world; // lambda2world[i] = world simulating set lambda i
int *world2lambda; // world2lambda[i] = lambda simulated by world i
int *world2root; // world2root[i] = root proc of world i
void print_status();
class FixGrem * fix_grem;
protected:
char *id_nh;
int pressflag;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Must have more than one processor partition to grem
Cannot use the grem command with only one processor partition. Use
the -partition command-line option.
E: Grem command before simulation box is defined
The grem command cannot be used before a read_data, read_restart, or
create_box command.
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: Tempering fix ID is not defined
The fix ID specified by the grem command does not exist.
E: Invalid frequency in grem command
Nevery must be > 0.
E: Non integer # of swaps in grem command
Swap frequency in grem command must evenly divide the total # of
timesteps.
E: Grem temperature fix is not valid
The fix specified by the grem command is not one that controls
temperature (nvt or npt).
E: Too many timesteps
The cummulative timesteps must fit in a 64-bit integer.
E: Grem could not find thermo_pe compute
This compute is created by the thermo command. It must have been
explicitly deleted by a uncompute command.
*/

View File

@ -28,9 +28,9 @@ class ComputePressure : public Compute {
public:
ComputePressure(class LAMMPS *, int, char **);
virtual ~ComputePressure();
void init();
double compute_scalar();
void compute_vector();
virtual void init();
virtual double compute_scalar();
virtual void compute_vector();
void reset_extra_compute_fix(const char *);
protected:

View File

@ -53,7 +53,7 @@ enum{ISO,ANISO,TRICLINIC};
FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg),
rfix(NULL), id_dilate(NULL), irregular(NULL), id_temp(NULL), id_press(NULL),
eta(NULL), eta_dot(NULL), eta_dotdot(NULL),
eta_mass(NULL), etap(NULL), etap_dot(NULL), etap_dotdot(NULL),
eta_mass(NULL), etap(NULL), etap_dot(NULL), etap_dotdot(NULL),
etap_mass(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal fix nvt/npt/nph command");
@ -442,7 +442,7 @@ etap_mass(NULL)
if (!atom->mu_flag)
error->all(FLERR,"Using update dipole flag requires atom attribute mu");
}
if ((tstat_flag && t_period <= 0.0) ||
(p_flag[0] && p_period[0] <= 0.0) ||
(p_flag[1] && p_period[1] <= 0.0) ||
@ -477,9 +477,9 @@ etap_mass(NULL)
// pre_exchange only required if flips can occur due to shape changes
if (flipflag && (p_flag[3] || p_flag[4] || p_flag[5]))
if (flipflag && (p_flag[3] || p_flag[4] || p_flag[5]))
pre_exchange_flag = 1;
if (flipflag && (domain->yz != 0.0 || domain->xz != 0.0 ||
if (flipflag && (domain->yz != 0.0 || domain->xz != 0.0 ||
domain->xy != 0.0))
pre_exchange_flag = 1;
}
@ -732,7 +732,7 @@ void FixNH::setup(int vflag)
t_current = temperature->compute_scalar();
tdof = temperature->dof;
// t_target is needed by NVT and NPT in compute_scalar()
// If no thermostat or using fix nphug,
// t_target must be defined by other means.
@ -1698,14 +1698,30 @@ void FixNH::reset_dt()
void *FixNH::extract(const char *str, int &dim)
{
dim=0;
if (strcmp(str,"t_target") == 0) {
if (tstat_flag && strcmp(str,"t_target") == 0) {
return &t_target;
} else if (strcmp(str,"mtchain") == 0) {
} else if (tstat_flag && strcmp(str,"t_start") == 0) {
return &t_start;
} else if (tstat_flag && strcmp(str,"t_stop") == 0) {
return &t_stop;
} else if (tstat_flag && strcmp(str,"mtchain") == 0) {
return &mtchain;
} else if (pstat_flag && strcmp(str,"mpchain") == 0) {
return &mtchain;
}
dim=1;
if (strcmp(str,"eta") == 0) {
if (tstat_flag && strcmp(str,"eta") == 0) {
return &eta;
} else if (pstat_flag && strcmp(str,"etap") == 0) {
return &eta;
} else if (pstat_flag && strcmp(str,"p_flag") == 0) {
return &p_flag;
} else if (pstat_flag && strcmp(str,"p_start") == 0) {
return &p_start;
} else if (pstat_flag && strcmp(str,"p_stop") == 0) {
return &p_stop;
} else if (pstat_flag && strcmp(str,"p_target") == 0) {
return &p_target;
}
return NULL;
}