Merge pull request #879 from mdelyser/USER-BOCS

USER-BOCS
This commit is contained in:
Steve Plimpton
2018-05-07 10:34:51 -06:00
committed by GitHub
17 changed files with 7321 additions and 2 deletions

View File

@ -101,7 +101,7 @@ set(DEFAULT_PACKAGES ASPHERE BODY CLASS2 COLLOID COMPRESS CORESHELL DIPOLE GRANU
KSPACE MANYBODY MC MEAM MISC MOLECULE PERI QEQ
REAX REPLICA RIGID SHOCK SNAP SRD)
set(OTHER_PACKAGES KIM PYTHON MSCG MPIIO VORONOI POEMS LATTE
USER-ATC USER-AWPMD USER-CGDNA USER-MESO
USER-ATC USER-AWPMD USER-BOCS USER-CGDNA USER-MESO
USER-CGSDK USER-COLVARS USER-DIFFRACTION USER-DPD USER-DRUDE USER-EFF
USER-FEP USER-H5MD USER-LB USER-MANIFOLD USER-MEAMC USER-MGPT USER-MISC
USER-MOFFF USER-MOLFILE USER-NETCDF USER-PHONON USER-QTB USER-REAXC USER-SMD

View File

@ -122,6 +122,7 @@ Package, Description, Doc page, Example, Library
Package, Description, Doc page, Example, Library
"USER-ATC"_#USER-ATC, atom-to-continuum coupling, "fix atc"_fix_atc.html, USER/atc, int
"USER-AWPMD"_#USER-AWPMD, wave-packet MD, "pair_style awpmd/cut"_pair_awpmd.html, USER/awpmd, int
"USER-BOCS"_#USER-BOCS, BOCS bottom up coarse graining, "fix bocs"_fix_bocs.html, USER/bocs, -
"USER-CGDNA"_#USER-CGDNA, coarse-grained DNA force fields, src/USER-CGDNA/README, USER/cgdna, -
"USER-CGSDK"_#USER-CGSDK, SDK coarse-graining model, "pair_style lj/sdk"_pair_sdk.html, USER/cgsdk, -
"USER-COLVARS"_#USER-COLVARS, collective variables library, "fix colvars"_fix_colvars.html, USER/colvars, int
@ -1625,6 +1626,43 @@ examples/USER/awpmd :ul
:line
USER-BOCS package :link(USER-BOCS),h4
[Contents:]
This package provides "fix bocs"_fix_bocs.html, a modified version
of "fix npt"_fix_nh.html which includes the pressure correction to
the barostat as outlined in:
N. J. H. Dunn and W. G. Noid, "Bottom-up coarse-grained models that
accurately describe the structure, pressure, and compressibility of
molecular liquids," J. Chem. Phys. 143, 243148 (2015).
[Authors:] Nicholas J. H. Dunn and Michael R. DeLyser (The Pennsylvania State University)
[Install or un-install:]
make yes-user-bocs
make machine :pre
make no-user-bocs
make machine :pre
[Supporting info:]
The USER-BOCS user package for LAMMPS is part of the BOCS software package:
"https://github.com/noid-group/BOCS"_https://github.com/noid-group/BOCS
See the following reference for information about the entire package:
Dunn, NJH; Lebold, KM; DeLyser, MR; Rudzinski, JF; Noid, WG.
"BOCS: Bottom-Up Open-Source Coarse-Graining Software."
J. Phys. Chem. B. 122, 13, 3363-3377 (2018).
Example inputs are in the examples/USER/bocs folder.
:line
USER-CGDNA package :link(USER-CGDNA),h4
[Contents:]

112
doc/src/fix_bocs.txt Normal file
View File

@ -0,0 +1,112 @@
<"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 bocs command :h3
[Syntax:]
fix ID group-ID bocs keyword values ... :pre
keyword = {temp} or {cgiso} or {analytic} or {linear_spline} or {cubic_spline}
{temp} values = Tstart Tstop Tdamp
{cgiso} values = Pstart Pstop Pdamp
{basis set}
{analytic} values = V_avg N_particles N_coeff Coeff_1 Coeff_2 ... Coeff_N
{linear_spline} values = input_filename
{cubic_spline} values = input_filename :pre
:ule
[Examples:]
fix 1 all bocs temp 300.0 300.0 100.0 cgiso 0.986 0.986 1000.0 analytic 66476.015 968 2 245030.10 8962.20 :pre
fix 1 all bocs temp 300.0 300.0 100.0 cgiso 0.986 0.986 1000.0 cubic_spline input_Fv.dat :pre
thermo_modify press 1_press :pre
[Description:]
These commands incorporate a pressure correction as described by
Dunn and Noid in "(Dunn1)"_#bocs-Dunn1 to the standard MTTK
barostat by Martyna et. al. in "(Martyna)"_#bocs-Martyna .
The first half of the command mimics a standard fix npt command:
fix 1 all bocs temp Tstart Tstop Tcoupl cgiso Pstart Pstop Pdamp :pre
The two differences are replacing {npt} with {bocs}, and replacing
{iso}/{aniso}/{etc} with {cgiso}.
The rest of the command details what form you would like to use for
the pressure correction equation. The choices are: {analytic}, {linear_spline},
or {cubic_spline}.
With either spline method, the only argument that needs to follow it
is the name of a file that contains the desired pressure correction
as a function of volume. The file should be formatted so each line has:
Volume_i, PressureCorrection_i :pre
Note both the COMMA and the SPACE separating the volume's
value and its corresponding pressure correction. The volumes in the file
should be uniformly spaced. Both the volumes and the pressure corrections
should be provided in the proper units, e.g. if you are using {units real},
the volumes should all be in cubic angstroms, and the pressure corrections
should all be in atomspheres. Furthermore, the table should start/end at a
volume considerably smaller/larger than you expect your system to sample
during the simulation. If the system ever reaches a volume outside of the
range provided, the simulation will stop.
With the {analytic} option, the arguments are as follows:
... analytic V_avg N_particles N_coeff Coeff_1 Coeff_2 ... Coeff_N :pre
Note that {V_avg} and {Coeff_i} should all be in the proper units, e.g. if you
are using {units real}, {V_avg} should be in cubic angstroms, and the
coefficients should all be in atmospheres * cubic angstroms.
[Restrictions:]
As this is computing a (modified) pressure, group-ID should be {all}.
The pressure correction has only been tested for use with an isotropic
pressure coupling in 3 dimensions.
By default, LAMMPS will still report the normal value for the pressure
if the pressure is printed via a {thermo} command, or if the pressures
are written to a file every so often. In order to have LAMMPS report the
modified pressure, you must include the {thermo_modify} command given in
the examples. For the last argument in the command, you should put
XXXX_press, where XXXX is the ID given to the fix bocs command (in the
example, the ID of the fix bocs command is 1 ).
This fix is part of the USER-BOCS 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:]
For more details about the pressure correction and the entire BOCS software
package, visit the "BOCS package on github"_bocsgithub and read the release
paper by Dunn et. al. "(Dunn2)"_#bocs-Dunn2 .
:link(bocsgithub,https://github.com/noid-group/BOCS)
:line
:link(bocs-Dunn1)
[(Dunn1)] Dunn and Noid, J Chem Phys, 143, 243148 (2015).
:link(bocs-Martyna)
[(Martyna)] Martyna, Tobias, and Klein, J Chem Phys, 101, 4177 (1994).
:link(bocs-Dunn2)
[(Dunn2)] Dunn, Lebold, DeLyser, Rudzinski, and Noid, J. Phys. Chem. B, 122, 3363 (2018).

View File

@ -20,6 +20,7 @@ Fixes :h1
fix_ave_time
fix_aveforce
fix_balance
fix_bocs
fix_bond_break
fix_bond_create
fix_bond_swap

View File

@ -135,6 +135,7 @@ fix_ave_histo.html
fix_ave_time.html
fix_aveforce.html
fix_balance.html
fix_bocs.html
fix_bond_break.html
fix_bond_create.html
fix_bond_react.html

View File

@ -0,0 +1,6 @@
This folder contains the files to run a NPT simulation of 1 site CG methanol
while employing a correction to the barostat.
The pair force was computed via the Multi-Scale Coarse-Graining method.
The resulting model was then iteratively pressure matched.
The model accurately reproduces both structural (RDF) and thermodynamic
(Pressure-Volume EoS) properties of the underlying OPLS-AA model of methanol.

View File

@ -0,0 +1,69 @@
units real
dimension 3
boundary p p p
atom_style atomic
newton on
timestep 1.0
read_data methanol.data
velocity all create 300.0 16802 dist gaussian
pair_style table spline 15000
pair_coeff 1 1 lammps_nb_MET-MET.table nb_METMET 12.0
neigh_modify delay 0 every 1 check yes one 10000
neighbor 12.0 bin
thermo 500
thermo_style custom step temp pe etotal press vol
variable STEP equal step
variable TEMP equal temp
## volume from cubic angstroms to cubic nm
variable VOL equal vol/1000.0
## pressure from atm to bar
variable PRESS equal press*1.01325
variable PXX equal pxx*1.01325
variable PYY equal pyy*1.01325
variable PZZ equal pzz*1.01325
variable PXY equal pxy*1.01325
variable PXZ equal pxz*1.01325
variable PYZ equal pyz*1.01325
## energy from kcal/mol to kJ/mol
variable KE equal ke*4.184
variable PE equal pe*4.184
variable UVDW equal evdwl*4.184
##### SPECIAL COMMANDS FOR FIX_BOCS #####
# ID group-ID style_name thermostat T_init T_end T_couple barostat P_start P_end P_couple pmatch_basis avg_vol N_sites N_coeffs coeff1 coeff2
fix 1 all bocs temp 300.0 300.0 100.0 cgiso 0.986 0.986 1000.0 analytic 66476.015 968 2 245030.10 8962.20
# Report the modified pressure
thermo_modify press 1_press
## Uncomment to save some data from simulation to files
#fix print_temp all print 500 "${STEP} ${TEMP}" file temp.dat screen no
#fix print_vol all print 500 "${STEP} ${VOL}" file vol.dat screen no
#fix print_press all print 500 "${STEP} ${PRESS}" file press.dat screen no
#fix print_ke all print 500 "${STEP} ${KE}" file kinetic_E.dat screen no
#fix print_pe all print 500 "${STEP} ${PE}" file potential_E.dat screen no
#fix print_ve all print 500 "${STEP} ${UVDW}" file vdw_E.dat screen no
#fix print_press_tens all print 500 "${STEP} ${PXX} ${PYY} ${PZZ} ${PXY} ${PXZ} ${PYZ}" file press_tens.dat screen no
#fix print_PV_eos all print 500 "${VOL} ${PRESS}" file pv_eos.dat screen no
## Prints a configuration to dump.txt every 500 steps
#dump 1 all custom 500 dump.txt id type x y z fx fy fz
# Write restart files to continue simulations
#restart 10000 state1.restart state2.restart
## Run for this many steps
run_style verlet
run 10000

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,143 @@
LAMMPS (20 Apr 2018)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:90)
using 1 OpenMP thread(s) per MPI task
units real
dimension 3
boundary p p p
atom_style atomic
newton on
timestep 1.0
read_data methanol.data
orthogonal box = (0 0 0) to (40.4635 40.4635 40.4635)
1 by 1 by 1 MPI processor grid
reading atoms ...
968 atoms
velocity all create 300.0 16802 dist gaussian
pair_style table spline 15000
pair_coeff 1 1 lammps_nb_MET-MET.table nb_METMET 12.0
WARNING: 78 of 2500 force values in table are inconsistent with -dE/dr.
Should only be flagged at inflection points (../pair_table.cpp:481)
neigh_modify delay 0 every 1 check yes one 10000
neighbor 12.0 bin
thermo 500
thermo_style custom step temp pe etotal press vol
variable STEP equal step
variable TEMP equal temp
## volume from cubic angstroms to cubic nm
variable VOL equal vol/1000.0
## pressure from atm to bar
variable PRESS equal press*1.01325
variable PXX equal pxx*1.01325
variable PYY equal pyy*1.01325
variable PZZ equal pzz*1.01325
variable PXY equal pxy*1.01325
variable PXZ equal pxz*1.01325
variable PYZ equal pyz*1.01325
## energy from kcal/mol to kJ/mol
variable KE equal ke*4.184
variable PE equal pe*4.184
variable UVDW equal evdwl*4.184
##### SPECIAL COMMANDS FOR FIX_BOCS #####
# ID group-ID style_name thermostat T_init T_end T_couple barostat P_start P_end P_couple pmatch_basis avg_vol N_sites N_coeffs coeff1 coeff2
fix 1 all bocs temp 300.0 300.0 100.0 cgiso 0.986 0.986 1000.0 analytic 66476.015 968 2 245030.10 8962.20
# Report the modified pressure
thermo_modify press 1_press
## Uncomment to save some data from simulation to files
#fix print_temp all print 500 "${STEP} ${TEMP}" file temp.dat screen no
#fix print_vol all print 500 "${STEP} ${VOL}" file vol.dat screen no
#fix print_press all print 500 "${STEP} ${PRESS}" file press.dat screen no
#fix print_ke all print 500 "${STEP} ${KE}" file kinetic_E.dat screen no
#fix print_pe all print 500 "${STEP} ${PE}" file potential_E.dat screen no
#fix print_ve all print 500 "${STEP} ${UVDW}" file vdw_E.dat screen no
#fix print_press_tens all print 500 "${STEP} ${PXX} ${PYY} ${PZZ} ${PXY} ${PXZ} ${PYZ}" file press_tens.dat screen no
#fix print_PV_eos all print 500 "${VOL} ${PRESS}" file pv_eos.dat screen no
## Prints a configuration to dump.txt every 500 steps
#dump 1 all custom 500 dump.txt id type x y z fx fy fz
# Write restart files to continue simulations
#restart 10000 state1.restart state2.restart
## Run for this many steps
run_style verlet
run 10000
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 10000, page size: 100000
master list distance cutoff = 24
ghost atom cutoff = 24
binsize = 12, bins = 4 4 4
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair table, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 4.691 | 4.691 | 4.691 Mbytes
Step Temp PotEng TotEng Press Volume
0 300 1061.5961 1926.3291 107.006 66250.679
500 314.54728 1034.1091 1940.7738 194.42689 65660.282
1000 301.41603 1030.7027 1899.5173 -91.966709 66262.543
1500 298.8308 1014.8276 1876.1905 -80.178606 67053.605
2000 294.78476 1046.8207 1896.521 50.592942 66316.735
2500 301.18564 1033.9214 1902.0719 40.48255 66607.667
3000 301.06632 1022.0381 1889.8447 47.582344 66341.947
3500 297.98361 989.80983 1848.7307 -204.69879 67462.078
4000 299.03493 1034.6571 1896.6083 89.188888 66457.385
4500 306.03351 985.4121 1867.5363 -51.102407 67519.446
5000 305.6903 1013.8613 1894.9963 -141.13704 67240.467
5500 292.23444 1029.5558 1871.905 20.764579 66683.876
6000 287.87735 1017.7325 1847.5226 -35.288049 66630.031
6500 305.26461 960.08118 1839.9891 -352.42596 67612.317
7000 300.34449 1055.0664 1920.7923 22.04027 66187.27
7500 305.48612 1038.6651 1919.2115 17.807254 66324.168
8000 316.03232 1034.6809 1945.6262 27.482857 66502.198
8500 294.28636 1038.8213 1887.085 -72.840559 66851.661
9000 316.69029 1065.7481 1978.5899 245.61677 65678.385
9500 297.46127 1034.5547 1891.97 54.23428 66892.627
10000 301.24799 1036.5432 1904.8735 7.7134029 66150.506
Loop time of 34.426 on 1 procs for 10000 steps with 968 atoms
Performance: 25.097 ns/day, 0.956 hours/ns, 290.478 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 33.324 | 33.324 | 33.324 | 0.0 | 96.80
Neigh | 0.12198 | 0.12198 | 0.12198 | 0.0 | 0.35
Comm | 0.42865 | 0.42865 | 0.42865 | 0.0 | 1.25
Output | 0.00059938 | 0.00059938 | 0.00059938 | 0.0 | 0.00
Modify | 0.42553 | 0.42553 | 0.42553 | 0.0 | 1.24
Other | | 0.1252 | | | 0.36
Nlocal: 968 ave 968 max 968 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 9112 ave 9112 max 9112 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 404392 ave 404392 max 404392 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 404392
Ave neighs/atom = 417.76
Neighbor list builds = 13
Dangerous builds = 0
Please see the log.cite file for references relevant to this simulation
Total wall time: 0:00:34

View File

@ -0,0 +1,143 @@
LAMMPS (20 Apr 2018)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:90)
using 1 OpenMP thread(s) per MPI task
units real
dimension 3
boundary p p p
atom_style atomic
newton on
timestep 1.0
read_data methanol.data
orthogonal box = (0 0 0) to (40.4635 40.4635 40.4635)
1 by 2 by 2 MPI processor grid
reading atoms ...
968 atoms
velocity all create 300.0 16802 dist gaussian
pair_style table spline 15000
pair_coeff 1 1 lammps_nb_MET-MET.table nb_METMET 12.0
WARNING: 78 of 2500 force values in table are inconsistent with -dE/dr.
Should only be flagged at inflection points (../pair_table.cpp:481)
neigh_modify delay 0 every 1 check yes one 10000
neighbor 12.0 bin
thermo 500
thermo_style custom step temp pe etotal press vol
variable STEP equal step
variable TEMP equal temp
## volume from cubic angstroms to cubic nm
variable VOL equal vol/1000.0
## pressure from atm to bar
variable PRESS equal press*1.01325
variable PXX equal pxx*1.01325
variable PYY equal pyy*1.01325
variable PZZ equal pzz*1.01325
variable PXY equal pxy*1.01325
variable PXZ equal pxz*1.01325
variable PYZ equal pyz*1.01325
## energy from kcal/mol to kJ/mol
variable KE equal ke*4.184
variable PE equal pe*4.184
variable UVDW equal evdwl*4.184
##### SPECIAL COMMANDS FOR FIX_BOCS #####
# ID group-ID style_name thermostat T_init T_end T_couple barostat P_start P_end P_couple pmatch_basis avg_vol N_sites N_coeffs coeff1 coeff2
fix 1 all bocs temp 300.0 300.0 100.0 cgiso 0.986 0.986 1000.0 analytic 66476.015 968 2 245030.10 8962.20
# Report the modified pressure
thermo_modify press 1_press
## Uncomment to save some data from simulation to files
#fix print_temp all print 500 "${STEP} ${TEMP}" file temp.dat screen no
#fix print_vol all print 500 "${STEP} ${VOL}" file vol.dat screen no
#fix print_press all print 500 "${STEP} ${PRESS}" file press.dat screen no
#fix print_ke all print 500 "${STEP} ${KE}" file kinetic_E.dat screen no
#fix print_pe all print 500 "${STEP} ${PE}" file potential_E.dat screen no
#fix print_ve all print 500 "${STEP} ${UVDW}" file vdw_E.dat screen no
#fix print_press_tens all print 500 "${STEP} ${PXX} ${PYY} ${PZZ} ${PXY} ${PXZ} ${PYZ}" file press_tens.dat screen no
#fix print_PV_eos all print 500 "${VOL} ${PRESS}" file pv_eos.dat screen no
## Prints a configuration to dump.txt every 500 steps
#dump 1 all custom 500 dump.txt id type x y z fx fy fz
# Write restart files to continue simulations
#restart 10000 state1.restart state2.restart
## Run for this many steps
run_style verlet
run 10000
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 10000, page size: 100000
master list distance cutoff = 24
ghost atom cutoff = 24
binsize = 12, bins = 4 4 4
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair table, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 3.37 | 3.37 | 3.371 Mbytes
Step Temp PotEng TotEng Press Volume
0 300 1061.5961 1926.3291 107.006 66250.679
500 314.54728 1034.1091 1940.7738 194.42689 65660.282
1000 301.41603 1030.7027 1899.5173 -91.966709 66262.543
1500 298.8308 1014.8276 1876.1905 -80.178606 67053.605
2000 294.78476 1046.8207 1896.521 50.592942 66316.735
2500 301.18564 1033.9214 1902.0719 40.482557 66607.667
3000 301.06631 1022.0381 1889.8447 47.582403 66341.947
3500 297.98353 989.81011 1848.7308 -204.69823 67462.076
4000 299.03465 1034.6603 1896.6108 89.196235 66457.338
4500 306.04532 985.37017 1867.5285 -51.094929 67519.735
5000 304.72903 1014.9543 1893.3184 -127.04402 67238.517
5500 292.52622 1025.6599 1868.8502 -19.753932 66716.551
6000 296.82719 1031.5184 1887.1059 -1.2609328 66368.611
6500 298.63312 1018.4299 1879.2229 -24.75835 66524.898
7000 303.25389 1005.9283 1880.0404 -96.273504 67349.674
7500 292.45089 1068.2863 1911.2595 103.23295 65778.08
8000 301.22765 1040.6294 1908.9011 -0.83635353 66831.038
8500 300.19765 1047.5856 1912.8883 -31.582343 66316.305
9000 295.1108 1023.8234 1874.4635 -88.165532 67192.344
9500 302.1087 1003.6348 1874.4459 -18.707065 66369.361
10000 296.3083 1004.126 1858.2178 -28.293045 66862.576
Loop time of 28.8053 on 4 procs for 10000 steps with 968 atoms
Performance: 29.994 ns/day, 0.800 hours/ns, 347.159 timesteps/s
95.2% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 19.929 | 21.765 | 23.391 | 27.8 | 75.56
Neigh | 0.067397 | 0.071231 | 0.077313 | 1.5 | 0.25
Comm | 3.9226 | 5.5183 | 7.3214 | 53.7 | 19.16
Output | 0.00069928 | 0.0016099 | 0.0043275 | 3.9 | 0.01
Modify | 1.0874 | 1.1376 | 1.1888 | 4.2 | 3.95
Other | | 0.3112 | | | 1.08
Nlocal: 242 ave 244 max 239 min
Histogram: 1 0 0 0 0 0 1 0 1 1
Nghost: 5718.5 ave 5736 max 5702 min
Histogram: 1 0 0 0 1 1 0 0 0 1
Neighs: 100703 ave 108064 max 93454 min
Histogram: 1 0 0 1 0 0 1 0 0 1
Total # of neighbors = 402813
Ave neighs/atom = 416.129
Neighbor list builds = 14
Dangerous builds = 0
Please see the log.cite file for references relevant to this simulation
Total wall time: 0:00:28

View File

@ -0,0 +1,991 @@
LAMMPS Description
968 atoms
0 bonds
0 angles
0 dihedrals
0 impropers
1 atom types
0 bond types
0 angle types
0 dihedral types
0 40.4635 xlo xhi
0 40.4635 ylo yhi
0 40.4635 zlo zhi
Masses
1 32.0424
Atoms
1 1 36.39 8.05 27.55
2 1 18.38 15.72 26.03
3 1 3.52 1.77 23.57
4 1 31.09 11.38 12.17
5 1 13.3 34.8 2.89
6 1 1.72 38.55 10.36
7 1 38.6 23.16 6.19
8 1 0.74 33.21 0.17
9 1 26.97 9.53 11.24
10 1 31.68 12.19 17.04
11 1 15.19 9.06 0.2
12 1 34.39 20.63 35.71
13 1 21.79 4.93 6.49
14 1 28.08 33.01 24.51
15 1 37.03 20.52 32.9
16 1 32.69 20 30.11
17 1 23.77 38.59 18.79
18 1 16.01 2.31 20.15
19 1 5.65 15.5 28.5
20 1 8.35 17.35 20.48
21 1 31.24 24.68 27.4
22 1 29.41 16.64 19.79
23 1 37.42 31.93 26.68
24 1 18.76 39.06 30.68
25 1 22.97 35 20.81
26 1 39.47 18.28 29.6
27 1 34.37 19.4 0.17
28 1 5.94 9.53 10.95
29 1 11.45 38.44 14.64
30 1 39.57 11.04 1.57
31 1 25.91 20.97 12.79
32 1 36.3 22.1 1.68
33 1 21.87 22.49 3.43
34 1 5.77 18.66 3.97
35 1 7 7.59 18.42
36 1 39.76 27.63 17.98
37 1 10.26 23.18 5.11
38 1 23.23 21.37 17.38
39 1 7.47 37.37 5.57
40 1 0.73 21.6 14.78
41 1 2.27 30.22 14.93
42 1 7.39 28.22 14.88
43 1 27.48 1.11 25.02
44 1 8.37 13.19 14.64
45 1 7.11 33.57 0.65
46 1 34.19 35.11 3.17
47 1 7.52 3.2 20.33
48 1 1.02 17.69 37.85
49 1 8.59 6.21 10.26
50 1 2.89 16.81 30.02
51 1 27.05 9.07 32.8
52 1 12.32 14.79 21.11
53 1 26.63 22.69 16.75
54 1 31.6 2.79 20.45
55 1 35.11 0.74 12.16
56 1 29.71 31.23 37.63
57 1 38.22 19.02 18.12
58 1 10.95 17.4 0.39
59 1 1.65 19.87 6.67
60 1 5.15 1.94 14.61
61 1 8.2 31.43 15.72
62 1 0.55 20.85 2.45
63 1 15.27 31.74 0.33
64 1 17.9 9.84 23.87
65 1 1.46 8.39 19.24
66 1 37.79 25.11 12.24
67 1 14.11 27.19 30.7
68 1 29.39 24.12 38.92
69 1 39.15 20.32 39.86
70 1 11.79 11.38 30.6
71 1 29.87 30.18 32.5
72 1 11.06 37.8 30.18
73 1 20.11 4.59 35.57
74 1 37.19 27.23 1.44
75 1 18.83 26.13 35.26
76 1 29.27 7.86 36.75
77 1 34.83 29.59 37.82
78 1 22.16 34.68 25.27
79 1 5.95 14.65 25.55
80 1 -0.04 3.48 29.91
81 1 14.69 6.75 38.06
82 1 16.8 23.69 3.14
83 1 14.52 38.72 21.61
84 1 37.36 15.13 27.98
85 1 29.79 17.03 27.55
86 1 21.01 0.39 0.9
87 1 35.3 6.45 14.41
88 1 13.74 34.69 10.01
89 1 13.83 7.92 31.88
90 1 26.88 37.81 22.36
91 1 31.76 7.34 10.03
92 1 11.98 37.49 19.61
93 1 29.19 31.74 20.69
94 1 39.15 10.97 21.06
95 1 1.47 5 22.21
96 1 -0.1 1.33 16.26
97 1 34.74 34.78 30.66
98 1 22.09 26.78 2.27
99 1 36.58 34.8 16.19
100 1 19.09 24.88 15.75
101 1 2.21 31.42 38.32
102 1 28.02 8.18 29.77
103 1 15.99 16.78 17.61
104 1 32.43 12.2 35.3
105 1 26.54 38.8 30.88
106 1 35.58 22.23 18.31
107 1 28.3 30.81 2.96
108 1 8.95 32.9 9.02
109 1 28.51 30.02 29.39
110 1 13 5.35 34.26
111 1 37.9 36.79 1.37
112 1 12.78 1.99 40.66
113 1 21.19 34.55 17.45
114 1 10.5 20 5.07
115 1 16.19 21.64 18.49
116 1 14.62 26.41 19.88
117 1 27.97 2.89 28.03
118 1 29.44 14.35 39.04
119 1 2.29 24.16 33.8
120 1 39.22 15.81 32.02
121 1 7.31 19.66 18.18
122 1 27.67 7.98 15.53
123 1 31.23 4.04 33
124 1 29.52 39.44 28.71
125 1 30.86 20.28 15.84
126 1 32.25 7.44 19.79
127 1 34.68 15.68 1.18
128 1 16.58 27.98 27.97
129 1 32.27 24.24 15.82
130 1 3.86 6.57 39.8
131 1 9.1 2.88 13.12
132 1 17.84 27.8 11.79
133 1 4.77 10.88 6.75
134 1 16.58 24.21 24.5
135 1 10.17 19.31 38.32
136 1 6.27 7.81 36.44
137 1 34.76 37.89 29.76
138 1 40.99 38.22 35.21
139 1 25.04 21.79 36.88
140 1 4.78 15.35 16.44
141 1 19.24 5.77 40.18
142 1 13.59 11.27 37.14
143 1 1.61 23.83 39.61
144 1 6.02 33.07 3.77
145 1 0.85 7.4 35.42
146 1 6.47 40.34 4.65
147 1 16.45 36.52 15.44
148 1 19.58 9.5 1.64
149 1 1.36 1.72 35.53
150 1 11.74 5.48 4.98
151 1 25.45 40.36 15.44
152 1 19.09 0.74 36.18
153 1 31.78 37.98 1.84
154 1 23.26 18.47 38.56
155 1 7.29 22.87 25.86
156 1 2.36 7.83 8.78
157 1 39.61 25.94 21.12
158 1 7.78 18.97 25.29
159 1 36.56 19.61 26.7
160 1 4.64 12.06 19.9
161 1 21.05 33.08 0.87
162 1 33.93 22.87 39.71
163 1 33.62 29.48 20.83
164 1 0.64 18.79 9.52
165 1 27.53 16.7 22.51
166 1 17.66 9.84 14.61
167 1 31.23 22.44 24.98
168 1 39.3 34.88 11.03
169 1 24.01 10.56 0.09
170 1 15.81 33.44 18.52
171 1 19.92 30.61 19.25
172 1 16.49 3.95 38.28
173 1 18.6 38.87 24.98
174 1 33.04 25.55 8.57
175 1 16.9 35.09 22.94
176 1 15.18 6.74 2.33
177 1 40.54 12.58 31.65
178 1 21.21 37.62 14.35
179 1 33.52 38.19 14.89
180 1 22.06 12.49 29.81
181 1 19.49 29.03 2.52
182 1 26.97 18.55 38.45
183 1 24.1 11.86 17.22
184 1 12.02 0.55 10.39
185 1 3.8 35.94 19.94
186 1 8.31 14.31 37.4
187 1 29.03 8.28 22.92
188 1 18.92 0.22 27.87
189 1 26.79 24.9 24.62
190 1 36.59 18.86 21.13
191 1 33.06 13.58 30.17
192 1 32.37 2.02 39.69
193 1 22.65 25.97 32.59
194 1 21.81 38.04 33.85
195 1 36.82 37.68 19.22
196 1 32.48 31.43 33.66
197 1 8.75 36.66 20.34
198 1 14.88 15.11 28.76
199 1 22.83 28.55 29.5
200 1 2.98 12.06 2.88
201 1 3.01 37 6.91
202 1 14.79 16.69 3.66
203 1 12.62 38.28 0.97
204 1 40.07 21.27 25.25
205 1 34.58 6.14 39.48
206 1 11.51 30.52 33.53
207 1 32.22 9.28 30.52
208 1 0.7 4.01 38.87
209 1 32.66 39.86 36.84
210 1 27.34 34.57 21.76
211 1 13.7 23.12 8.34
212 1 28 9.71 6.29
213 1 31.69 26.39 38.19
214 1 2.04 20.93 36.95
215 1 30.31 26.61 12.31
216 1 25.37 20.63 9.4
217 1 16.29 5.62 14.04
218 1 1.36 1.95 4.97
219 1 38.75 26.8 5.62
220 1 40.28 29.89 12.27
221 1 0.52 31.22 7.58
222 1 15.64 21.64 14.71
223 1 10.11 28.8 9.67
224 1 19.03 37.3 38.87
225 1 12.07 6.46 19.25
226 1 36.22 21.11 15.02
227 1 12.49 15.69 7.21
228 1 25.2 30.8 19.71
229 1 8.24 35.22 23.02
230 1 9.94 4.1 1.86
231 1 4.31 19.07 26.23
232 1 22.03 19.99 14.19
233 1 30.84 5.24 2.75
234 1 35.51 30.52 32.3
235 1 3.88 20.41 39.48
236 1 30.86 0.39 31.56
237 1 25.91 27.12 39.04
238 1 33.39 5.37 4.51
239 1 20.58 28.33 37.53
240 1 11.83 21.81 39.11
241 1 39.94 5.2 5.73
242 1 29.89 33.21 15.44
243 1 23.42 15.38 30.98
244 1 28.66 11.56 16.13
245 1 16.74 20.43 21.72
246 1 9.31 0.26 22.42
247 1 4.38 37.78 14.12
248 1 13.86 4.65 16.99
249 1 6.72 27.43 8.14
250 1 33.45 18.74 19.98
251 1 17.34 14.55 35.57
252 1 14.21 37.39 28.76
253 1 26.29 29.06 22.26
254 1 36.52 26.18 15.62
255 1 17.6 8.81 10.14
256 1 21.95 25.19 17.69
257 1 17.4 7.46 18.59
258 1 6.15 29.7 35.99
259 1 9.73 37.39 24.96
260 1 28.58 28.21 36.19
261 1 8.24 13.72 32.2
262 1 13.25 34.6 22.71
263 1 38.58 7.98 10.3
264 1 18.29 29.75 23.51
265 1 3.74 13.22 31.05
266 1 30.78 0.7 24.57
267 1 24.5 12.94 10.81
268 1 33.88 26.44 1.96
269 1 29.92 22.6 9.59
270 1 15.87 21.94 33.44
271 1 23.63 8.41 10.95
272 1 26.53 0.01 7.58
273 1 33.88 8.73 2.52
274 1 39.74 31.46 36.87
275 1 13.28 39.54 37
276 1 24.74 35.41 32.14
277 1 17.3 34.32 6.85
278 1 24.43 13.36 5.63
279 1 31.23 17.12 7.11
280 1 15.58 37.85 34.06
281 1 25.93 38.17 3.13
282 1 31.84 34.34 39.67
283 1 14.41 17.74 14.43
284 1 9.4 1.78 29.53
285 1 8.29 29.61 6.81
286 1 7.13 17.36 -0.19
287 1 38.04 19.52 11.25
288 1 34.87 32.77 38.75
289 1 27.18 22.49 5.69
290 1 29.69 27.11 20.29
291 1 25.23 2.22 39.75
292 1 35.8 17.1 36.09
293 1 20.53 17.67 35.47
294 1 23.41 31.36 2.24
295 1 25.17 20.09 29.2
296 1 0.32 9.28 39.74
297 1 22.49 18.77 29.89
298 1 38.64 0.85 39.43
299 1 18.55 5.87 9.15
300 1 35.21 28.69 25.76
301 1 39.24 15.68 6.36
302 1 5.91 13.29 3.92
303 1 0.5 21.25 21.62
304 1 34.39 9.43 18.59
305 1 22.7 1.03 8.79
306 1 40.15 16.25 24.87
307 1 16.52 13.69 10.25
308 1 6.88 36.77 39.57
309 1 11.95 37.1 7.59
310 1 19.49 26.66 9.76
311 1 36.34 29.89 29.31
312 1 18.76 35.75 12.72
313 1 23.15 38.68 4.96
314 1 10.78 17.11 34.38
315 1 13.58 14.89 0.34
316 1 4.73 33.88 32.41
317 1 12.68 29.51 17.07
318 1 9.45 27 4.98
319 1 29.02 19.27 6.45
320 1 19.09 36.35 19.84
321 1 23.13 6.75 19.78
322 1 8.02 34.1 26.95
323 1 26.03 14.18 1.57
324 1 15.01 34.16 28.32
325 1 28.71 36.53 15.08
326 1 9.64 13.94 24.05
327 1 1.24 36.96 30.41
328 1 15.44 31.39 4.77
329 1 25.72 14.42 38.01
330 1 7.35 31.84 12.61
331 1 32.33 16.79 3.07
332 1 15.86 0.39 8.39
333 1 27.69 7.32 25.74
334 1 27.77 1.17 4.96
335 1 29.38 26.68 23.48
336 1 14.95 3.41 2.13
337 1 2.45 38.8 17.54
338 1 9.18 27.58 31.12
339 1 36.62 15.23 23.66
340 1 12.23 30.74 10.85
341 1 12.87 0.23 24.05
342 1 33.85 35.75 19.95
343 1 36.36 30.35 13.94
344 1 38.11 8.15 5.49
345 1 15.58 1.29 29.89
346 1 25.95 30.47 39.19
347 1 27.69 15.5 13.37
348 1 25.48 13.7 25.29
349 1 0.72 39.53 4.41
350 1 8.75 15.21 5.79
351 1 10.49 26.46 27.25
352 1 16.9 20.08 26.58
353 1 3.95 4.33 20.25
354 1 18.03 7.66 31.87
355 1 21.35 2.88 17.26
356 1 32.92 22.68 6.56
357 1 21.72 4.62 12.02
358 1 37.6 6.04 25.2
359 1 22.22 23.78 7.36
360 1 12.53 19.53 10.6
361 1 17.87 26.05 19.32
362 1 20.94 30.75 31.23
363 1 22.33 11.87 35.79
364 1 28.54 6.88 10.52
365 1 29.58 26.13 8.56
366 1 19.06 24.2 21.73
367 1 25.05 38.41 26.91
368 1 18.66 40.35 4.69
369 1 11.87 6.75 12.72
370 1 17.99 11.82 37.79
371 1 8.97 24.7 0.61
372 1 26.89 1.52 31.33
373 1 6.88 22.51 14.06
374 1 29.51 3.87 13.35
375 1 0.95 14.53 10
376 1 26.43 31.54 11.63
377 1 6.1 35.16 11.28
378 1 10.2 15.91 28.7
379 1 13.24 25.67 24.3
380 1 34.5 3.91 20.22
381 1 22.48 2.17 5.2
382 1 24.26 16.15 11.3
383 1 20.63 32.57 27.3
384 1 39.95 6.81 1.05
385 1 24.38 34.78 38.27
386 1 4.63 23.56 0.55
387 1 33.5 8.38 36.97
388 1 5.76 27.14 12.06
389 1 37.45 26.58 32.25
390 1 2.8 8.96 32.35
391 1 5.29 39.84 30.1
392 1 29.2 26.49 31.18
393 1 33.91 27.39 12.89
394 1 3.37 14.19 39.34
395 1 30.68 28.7 27
396 1 2.59 18.45 19.38
397 1 13.55 0.75 27.35
398 1 3.82 20.12 16.1
399 1 37.72 0.72 34.29
400 1 23.85 3.32 19.46
401 1 4.78 0.67 38.23
402 1 22.78 23.26 38.6
403 1 11.56 39.29 4.23
404 1 21.38 32.45 5.39
405 1 8.32 19.42 9.52
406 1 28.43 31.07 17.8
407 1 11.02 5.75 39.84
408 1 27.36 36.71 7.62
409 1 34.22 16.74 27.88
410 1 3.22 22.01 27.42
411 1 29.2 15.76 32.33
412 1 25.29 23.44 2.2
413 1 10.8 32.43 39.96
414 1 32.2 1.41 4.44
415 1 32.94 15.59 37.21
416 1 6.8 8.08 0.84
417 1 10.42 9.91 37.73
418 1 1.18 31.39 3.93
419 1 10.1 36.38 38.57
420 1 32.89 26.96 35.07
421 1 28.12 11.93 25.96
422 1 4.9 29.15 -0.25
423 1 2.21 27.99 3.72
424 1 11.33 3.94 25.55
425 1 3.3 30.26 10.78
426 1 11.57 27.26 19.31
427 1 21.79 32.67 13.3
428 1 4.96 26.53 33.78
429 1 33.41 32.87 18.46
430 1 13.92 30.37 20.3
431 1 16.91 3.5 11.57
432 1 -0.06 4.42 34.05
433 1 7.04 24.33 16.85
434 1 28.66 11.93 19.73
435 1 30.21 1.75 36.36
436 1 3.91 6.2 6.26
437 1 7.01 25.6 27.34
438 1 34 1.06 18.12
439 1 29.14 8.5 3.09
440 1 40.13 23.52 16.94
441 1 21.69 22.38 26.78
442 1 18.44 32.85 39.84
443 1 38.87 1.99 27.63
444 1 10.47 12.02 33.68
445 1 9.65 19.94 21.88
446 1 25.04 8.12 27.89
447 1 12.18 16.78 31.09
448 1 38.31 8.2 30.6
449 1 0.11 5.21 18.12
450 1 23.48 7.88 38.7
451 1 7.9 11.99 7.04
452 1 2.09 34.5 8.64
453 1 19 18.3 0.42
454 1 37.94 18.06 15.21
455 1 25.05 33.33 29.55
456 1 33.74 37.88 34.54
457 1 36.07 29.5 17.47
458 1 27.7 13.54 10.46
459 1 27.18 23.18 21.06
460 1 17.49 2.04 23.22
461 1 23.84 26.31 36.76
462 1 31.08 32.24 24.92
463 1 17.7 28.73 16.66
464 1 1.59 15.5 18.04
465 1 33.49 18.99 8.75
466 1 1.19 8.43 28.01
467 1 5.43 20.26 22.25
468 1 23.3 5.31 15.05
469 1 27.09 4.47 21.96
470 1 26.61 27.48 28.83
471 1 13 3.18 13.4
472 1 36.09 34.56 10.58
473 1 19.09 15.33 7.48
474 1 12.94 18.16 22.07
475 1 38.8 16.16 36.73
476 1 16.73 39.94 -0.46
477 1 3.64 12.56 9.42
478 1 27.64 3.31 18.29
479 1 15.58 12.36 4.52
480 1 15.37 11.18 18.19
481 1 37.52 11.23 14.94
482 1 37.16 2.3 10.01
483 1 36.19 9.86 21.99
484 1 36.38 21.09 4.76
485 1 38.15 12.94 25.35
486 1 15.1 29.65 24.64
487 1 16.54 29.12 38.85
488 1 33.86 11.39 14.4
489 1 27.99 18.55 9.87
490 1 0.66 1.26 8.28
491 1 29.89 33.84 29.75
492 1 6.64 33.3 6.74
493 1 31.47 0.39 11.25
494 1 8.76 15.96 10.74
495 1 39.6 15.8 21.04
496 1 22.38 28.55 19.85
497 1 25.87 6.23 5.4
498 1 36.24 26.67 38.61
499 1 23.05 8.79 7.01
500 1 32.05 4.16 8.87
501 1 35.3 13.36 38.03
502 1 39.91 25.91 36.45
503 1 32.17 17.27 31.82
504 1 29.99 21.54 20.57
505 1 9.39 0.57 34.17
506 1 22.12 17.22 13.03
507 1 15.23 16.27 24.39
508 1 26.32 25.89 13.26
509 1 39.25 3.46 1.34
510 1 32.56 10.28 7.96
511 1 25.76 14.8 34.79
512 1 32.12 5.38 36.94
513 1 17.74 15.37 14.87
514 1 21.93 3.26 25.26
515 1 24.05 0.48 36.27
516 1 8.2 19.02 33.92
517 1 33.07 25.88 24.79
518 1 12.54 0.52 32.8
519 1 18.5 6.34 23.21
520 1 35.93 10.26 34.79
521 1 19.33 11.86 6.94
522 1 1.63 5.31 25.43
523 1 30.62 36.78 30.77
524 1 12.25 26.79 13.74
525 1 21 1.47 13.99
526 1 22.24 29.36 34.11
527 1 29.05 1.74 0.24
528 1 19.34 25.14 39.77
529 1 12.9 25.42 35.42
530 1 4.1 20.69 33.31
531 1 35.38 5.64 9.17
532 1 5.01 34.03 23.36
533 1 15.7 10.01 7.1
534 1 25.95 19.48 24.73
535 1 11.15 4.16 21.43
536 1 1.47 3.23 13.72
537 1 26.54 36.76 17.48
538 1 7.65 30.6 27.73
539 1 18.59 3.05 3.32
540 1 6.92 36.14 15.67
541 1 12.72 36.21 32.72
542 1 16.33 26.24 14.35
543 1 21.66 12.61 26.18
544 1 15.72 32.37 10.71
545 1 27.57 15.75 4.97
546 1 20.18 39.13 18.44
547 1 18.77 17.09 4.74
548 1 12.75 14.08 17.89
549 1 10.71 37.26 35.13
550 1 20.88 32.12 9.94
551 1 8.69 0.06 9.54
552 1 1.58 12.61 26.42
553 1 10.29 22.73 13.14
554 1 23.31 25.7 4.97
555 1 36.7 30.78 5.26
556 1 5.6 28.71 23.76
557 1 13.04 7.91 15.46
558 1 40.24 33.53 34.5
559 1 39.7 18.36 34.25
560 1 23.52 22.84 21.14
561 1 32.82 27.9 5.02
562 1 1.52 33.77 27.04
563 1 31.48 18.68 13.1
564 1 34 24.07 34.55
565 1 34.72 25.03 30.45
566 1 18.92 5.7 16.54
567 1 15.26 40.2 16.97
568 1 24.55 16.8 26.26
569 1 17.59 18.03 31.75
570 1 12.08 27.38 38.42
571 1 31.42 21.91 35.27
572 1 9.68 6.93 36.96
573 1 22.93 15.19 2.53
574 1 11.83 24.38 2.24
575 1 16.02 12.44 31.42
576 1 20.24 10.7 20.15
577 1 37.38 17.37 40.63
578 1 3.74 17.41 22.49
579 1 23.21 9.86 23.99
580 1 2.49 36.56 38
581 1 20.51 1.54 22.62
582 1 8.23 33.44 33.15
583 1 30.93 8.12 14.85
584 1 36.35 3.92 37.65
585 1 5.57 10.63 28.64
586 1 30.24 22.84 4.6
587 1 3.24 15.97 7.43
588 1 22.64 0.7 27.09
589 1 1.64 40.37 27.74
590 1 30.53 35.28 3.5
591 1 25.55 4.7 2.53
592 1 22.04 29.56 12.29
593 1 3.42 8.76 2.28
594 1 16.37 37.65 7.87
595 1 20.65 12.75 2.9
596 1 8.07 20.49 2.6
597 1 20.04 24.72 28.2
598 1 29.75 0.33 16.1
599 1 31.36 33.09 8.02
600 1 8.6 10.56 2.65
601 1 20.55 33.93 33.71
602 1 23.11 -0.42 31.7
603 1 27.07 6.69 18.5
604 1 15.65 23.71 11.8
605 1 23.88 24.17 30.1
606 1 0.79 32.46 16.99
607 1 25.96 20.19 2.43
608 1 35.16 12.33 4.73
609 1 16.53 17.12 10.45
610 1 19.74 35.94 8.62
611 1 -0.24 28.29 28.3
612 1 33.15 4.39 12.64
613 1 24.27 18.59 33.59
614 1 0.3 9.61 23.54
615 1 5.52 10.89 34.71
616 1 9.85 26.26 34.42
617 1 25.54 37.74 39.13
618 1 4.11 29.65 5.81
619 1 30.86 22.72 31.07
620 1 39.54 34.54 23.24
621 1 15.03 9.34 12.11
622 1 21.93 38.49 9.53
623 1 18.21 19.2 16.65
624 1 38.71 32.83 14.49
625 1 37.74 36.38 31.43
626 1 7.57 2.2 1.07
627 1 12.27 10.05 1.53
628 1 18.32 19.55 12.51
629 1 12.49 24.81 16.57
630 1 35.07 31.62 22.85
631 1 4.91 19.68 30.45
632 1 11.36 2.92 36.19
633 1 26.76 35.53 1.29
634 1 1.13 14.55 14.77
635 1 8.2 25.14 37.83
636 1 -0.5 0.85 23.08
637 1 17.95 18.01 37.61
638 1 2.44 37.75 0.71
639 1 3.77 23.41 17.91
640 1 3.34 14.41 35.96
641 1 14.89 31.42 31.76
642 1 15.72 22.91 30.66
643 1 37.42 3.13 5.51
644 1 1.51 38.14 20.78
645 1 12.55 18.16 25.81
646 1 18.98 15.88 29.21
647 1 17.33 32.04 27.87
648 1 31.21 37.38 22.84
649 1 11.88 32.4 31
650 1 25.1 2.19 12.69
651 1 38.85 39.1 27.87
652 1 5.9 16.14 33.34
653 1 6.3 17.5 36.46
654 1 2.28 25.02 12.53
655 1 19.33 25.49 6.78
656 1 8.57 37.22 8.4
657 1 22 38.61 24.14
658 1 35.63 3.23 2.15
659 1 37.35 22.64 9.99
660 1 25.46 28.22 33.18
661 1 13.99 4.21 9.07
662 1 9.04 40.64 37.59
663 1 17.46 6.49 27.93
664 1 1.42 6.14 30.2
665 1 25.82 28.28 10.33
666 1 30.71 30.2 -0.02
667 1 8.05 29.32 0.62
668 1 6.49 3.17 6.65
669 1 7.77 39.88 26.03
670 1 28.53 12.14 30.96
671 1 6.25 12.64 39.55
672 1 3.78 7.96 25.83
673 1 28.69 29.58 6.61
674 1 5.32 18.55 8.6
675 1 5.3 11.44 14.19
676 1 39 12.03 36.99
677 1 11.08 32.51 19.12
678 1 38.48 1.03 20.34
679 1 38.47 11.31 10.12
680 1 3.65 20.3 11.99
681 1 27.18 11.39 3.4
682 1 29.61 20.72 38.86
683 1 31.88 4.82 24.75
684 1 14.52 36.61 25.66
685 1 21.68 19.96 24.58
686 1 1.2 13.59 23.16
687 1 37.71 5.97 35.28
688 1 19.81 22.04 37.54
689 1 40.67 27.07 10.14
690 1 38.57 13.96 13.75
691 1 18.15 17.19 21.63
692 1 3.76 27.51 30.32
693 1 11.39 25.09 7.96
694 1 28.95 22.02 13.33
695 1 31.72 35.46 35.33
696 1 15.17 6.03 23.94
697 1 9.63 9.42 22.18
698 1 27.76 4.1 38.61
699 1 27.37 37.46 11.07
700 1 35.69 15.24 6.39
701 1 21.92 33.09 36.84
702 1 4.79 15.51 13.28
703 1 9.31 26.44 22.8
704 1 30.14 38.31 6.98
705 1 4.33 6.66 13.77
706 1 3.66 27.15 37.83
707 1 0.44 24.06 29.62
708 1 35.4 9.59 9.8
709 1 33.77 39.44 21.41
710 1 29.2 18.74 24.38
711 1 34.14 30 7.36
712 1 16.02 4.37 6.43
713 1 35.75 39.41 40.58
714 1 15.32 11.84 23.24
715 1 32.52 24.39 20.97
716 1 11.45 2.9 7.16
717 1 7.07 34.28 36.04
718 1 11.41 8.83 7.57
719 1 29.29 27.79 16.55
720 1 2.41 9.16 16.14
721 1 13.14 38.07 12.01
722 1 37.45 14.26 18.3
723 1 34.64 18.37 16.65
724 1 29.26 34.28 33.46
725 1 13.67 19.83 3.04
726 1 19.64 2.6 32.95
727 1 23.61 23.55 11.34
728 1 6.87 22.5 36.5
729 1 36.78 28.34 22.96
730 1 11.3 13.98 9.84
731 1 20.82 19.73 9.16
732 1 12.65 12.61 6.16
733 1 9.25 7.81 25.28
734 1 25.4 29.36 5.84
735 1 17.86 9.06 4.45
736 1 35.84 40.51 26.06
737 1 26.97 21.28 31.59
738 1 39.94 38.14 24.86
739 1 3.16 17.07 2.49
740 1 14.08 15.4 36.99
741 1 26.1 13.88 29.68
742 1 26.36 20 20.56
743 1 12.02 6.35 27.65
744 1 11.21 19.14 15.27
745 1 35.86 22.93 26.81
746 1 32.26 12.62 2.68
747 1 29.36 4.91 6.86
748 1 20.14 7.97 29.32
749 1 25.53 11.2 36.18
750 1 30.36 14.25 24.79
751 1 29.84 38.37 39.29
752 1 15.79 36.4 4.26
753 1 32.59 14.17 10.29
754 1 13.89 34.54 15.42
755 1 12.13 33.62 7.27
756 1 25.86 23.81 33.97
757 1 18.4 31.27 34.96
758 1 6.58 40.49 17.57
759 1 5.64 39.86 23.05
760 1 25.32 33 16.24
761 1 0.98 10.99 12.21
762 1 32.86 23.75 12.41
763 1 32.91 1.3 26.94
764 1 9.46 8.66 31.47
765 1 17.49 15.53 1.82
766 1 8.17 18.16 15.3
767 1 4.84 30.63 26.32
768 1 6.75 37.1 30.8
769 1 7.16 5.72 15.74
770 1 20.09 17.82 19.18
771 1 1.02 27.99 32.54
772 1 21.44 1.28 38.4
773 1 21.06 14.62 37.08
774 1 27.82 18.99 15.96
775 1 33.51 21.46 22.82
776 1 8.29 2.79 17.09
777 1 18.14 11.41 28.62
778 1 17.94 28.54 32.73
779 1 36.1 9.54 40.1
780 1 36.55 2.62 22.97
781 1 27.29 10.41 39.35
782 1 22.04 37.76 2.02
783 1 23.01 30.08 16.58
784 1 6.34 15.31 21.61
785 1 7.3 22.39 7.74
786 1 38 5.77 21.04
787 1 32.93 12.85 26.58
788 1 22.43 15.36 16.72
789 1 18.54 2.38 8.35
790 1 0.62 31.88 23.84
791 1 39.56 30.3 21.7
792 1 7.56 11.78 25.2
793 1 16.05 3.98 26.72
794 1 24.33 36.16 13.16
795 1 26.38 31.24 35.77
796 1 14.94 26.46 -0.44
797 1 3.53 28 19.79
798 1 11.94 32.27 25.88
799 1 34.59 15.06 19.23
800 1 3.51 10.62 24.22
801 1 13.18 21.59 35.52
802 1 19.95 13.15 16.23
803 1 24.48 9.63 19.58
804 1 37.7 34.42 19.19
805 1 22.84 19.59 5.52
806 1 34.62 36.75 12.26
807 1 10.33 21.85 25.15
808 1 21.51 29.7 7.23
809 1 34.05 3 30.58
810 1 12.43 29.67 5.59
811 1 18.94 21.82 29.31
812 1 39.21 17.7 4
813 1 7.04 2.97 37.09
814 1 30.87 10.8 22.09
815 1 33.54 7.16 23.15
816 1 32.82 32.27 11.47
817 1 21.8 14.54 21.61
818 1 38.83 6.69 13.68
819 1 30.37 39.05 18.66
820 1 31.21 15.08 13.79
821 1 1.07 11.97 35.04
822 1 26.02 35.8 35.47
823 1 27.64 33.32 8.04
824 1 39.57 36.88 4.39
825 1 14.48 30.97 36.87
826 1 36.86 27.04 8.96
827 1 34.73 26.27 19.55
828 1 29.24 29.12 10.86
829 1 8.35 39.9 13.92
830 1 5.43 5.1 31.85
831 1 18.36 22.31 5.89
832 1 2.46 27.35 16.36
833 1 21.76 37.05 27.45
834 1 29.32 17.11 2.28
835 1 3.58 21.98 3.09
836 1 10.61 9.17 10.69
837 1 25.15 17.54 6.05
838 1 3.4 22.81 7.03
839 1 16.12 23.66 38.29
840 1 8.67 4.79 31.8
841 1 27.82 0.02 20.45
842 1 4.46 2.44 29.46
843 1 39.11 22.69 36.33
844 1 4.26 2.45 34.43
845 1 13.47 3.5 30.9
846 1 26.92 22.55 26.54
847 1 5.66 30.07 30.65
848 1 31.58 35.68 25.46
849 1 31.23 5.65 27.78
850 1 8.76 19.49 28.61
851 1 24.13 8.02 2.3
852 1 8.27 6.8 4.93
853 1 15.31 34.36 37.63
854 1 3.96 35.82 26.97
855 1 24.95 4.01 29.96
856 1 26.54 25.95 19.15
857 1 17.76 29.23 5.35
858 1 16.95 21.79 0.78
859 1 17.52 39.02 13.01
860 1 20.5 22.86 10.41
861 1 27.67 27.07 2.78
862 1 37.1 31.02 2.03
863 1 37.17 2.24 16.34
864 1 24.51 34.69 10.28
865 1 17.38 13.41 21.4
866 1 38 9.37 17.94
867 1 35.67 25.27 5.3
868 1 22.66 4.47 0.97
869 1 20.06 28.17 27.42
870 1 27.33 2.38 9.79
871 1 36.4 10.73 31.42
872 1 0.53 9.45 6.96
873 1 12.89 9.03 25.13
874 1 37.5 23.54 21.33
875 1 7.19 28.87 19.11
876 1 21.45 32.81 22.71
877 1 8.76 8.3 13.51
878 1 10.24 31.46 3.73
879 1 15.93 9.7 35.82
880 1 14.92 19.46 39.91
881 1 13.02 22.03 27.3
882 1 34.57 1.59 34.73
883 1 0.5 21.04 30.77
884 1 24.84 32.67 4.97
885 1 30.96 4.24 17.34
886 1 23.63 8.84 32.7
887 1 6.84 3.83 26.13
888 1 31.57 9.23 27.06
889 1 14.99 35.15 0.07
890 1 36.29 40.01 8.22
891 1 18 36.86 35.66
892 1 11.5 29.22 0.37
893 1 35.18 36.03 24.88
894 1 14.49 9.41 27.73
895 1 10.12 12.91 0.99
896 1 11.91 29.4 27.3
897 1 25.51 28.68 14.56
898 1 29.32 12.95 35.59
899 1 2.96 0.59 0.53
900 1 5.47 5.35 9.88
901 1 8.93 39.33 1.78
902 1 39.68 38.49 12.33
903 1 3.78 1.24 11.85
904 1 29.54 33.76 11.99
905 1 7.48 6.16 28.8
906 1 10.83 33.63 36.96
907 1 22.03 19.06 21.39
908 1 31.68 31 3.71
909 1 35.32 0.3 30.88
910 1 9.93 35.42 14.85
911 1 21.01 10.66 12.97
912 1 22.54 21.62 33.81
913 1 12.48 12.46 25.15
914 1 0.64 27.93 39.01
915 1 13.88 27.29 3.8
916 1 13.03 20.93 19.25
917 1 5.7 37.48 34.09
918 1 2.27 35.08 14.29
919 1 21.75 26.44 13.48
920 1 10.66 11.35 17.49
921 1 36.67 34.72 27.69
922 1 8.23 22.8 33.36
923 1 4.79 4.05 3.48
924 1 35.61 35 35.21
925 1 31.38 37.67 10.06
926 1 26.44 35.58 27.52
927 1 25.53 4.47 8.03
928 1 38.81 33.08 31.55
929 1 29.62 39.68 34.32
930 1 0.76 3.88 10.61
931 1 24.06 30.13 24.84
932 1 18.77 8.38 37.02
933 1 2.81 24.14 24.45
934 1 11.33 1.25 19.69
935 1 1.61 26.53 6.48
936 1 9.52 29.81 35.98
937 1 7.83 31.75 20.77
938 1 9.67 23.93 20.57
939 1 33.96 16.95 24.79
940 1 37.55 37.31 37.17
941 1 6.2 24.19 31.13
942 1 17.79 4.35 30.56
943 1 13.51 20.25 30.45
944 1 14.76 12.48 13.13
945 1 32.56 27.52 30.41
946 1 35.81 40.1 4.85
947 1 39.01 25.43 27.07
948 1 22.8 5.53 32.22
949 1 2.35 39.55 32.5
950 1 3.31 31.51 34.19
951 1 40.63 37.01 15.9
952 1 18.43 35.86 30.77
953 1 37.01 39.18 15.57
954 1 6.74 6.26 21.64
955 1 5.22 24.94 22.53
956 1 25.2 11.08 13.98
957 1 14.66 -0.14 5
958 1 14.96 8.38 20.32
959 1 1.89 23.03 10.1
960 1 33.01 11.23 39.32
961 1 32.68 6.79 32.04
962 1 6.25 25.81 3.82
963 1 18.37 31.97 14.57
964 1 9.74 30.21 22.71
965 1 3.09 1.52 19.38
966 1 33.98 37.34 5.63
967 1 32.19 28.35 15.8
968 1 24.4 15.49 19.39

View File

@ -57,7 +57,7 @@ PACKAGE = asphere body class2 colloid compress coreshell dipole gpu \
molecule mpiio mscg opt peri poems \
python qeq reax replica rigid shock snap srd voronoi
PACKUSER = user-atc user-awpmd user-cgdna user-cgsdk user-colvars \
PACKUSER = user-atc user-awpmd user-bocs user-cgdna user-cgsdk user-colvars \
user-diffraction user-dpd user-drude user-eff user-fep user-h5md \
user-intel user-lb user-manifold user-meamc user-meso \
user-mgpt user-misc user-molfile \

19
src/USER-BOCS/README Normal file
View File

@ -0,0 +1,19 @@
This user package implements the pressure correction to the barostat as
outlined in:
N. J. H. Dunn and W. G. Noid, "Bottom-up coarse-grained models that
accurately describe the structure, pressure, and compressibility of
molecular liquids," J. Chem. Phys. 143, 243148 (2015).
doi: 10.1063/1.4937383
The USER-BOCS user package for LAMMPS is part of the BOCS software package:
https://github.com/noid-group/BOCS
See the following reference for information about the entire package:
Dunn, NJH; Lebold, KM; DeLyser, MR; Rudzinski, JF; Noid, WG.
"BOCS: Bottom-Up Open-Source Coarse-Graining Software."
J. Phys. Chem. B. 122, 13, 3363-3377 (2018).
Example inputs are in the examples/USER/bocs folder.

View File

@ -0,0 +1,444 @@
/* ----------------------------------------------------------------------
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.
-------------------------------------------------------------------------
USER-BOCS written by: Nicholas J. H. Dunn and Michael R. DeLyser
from The Pennsylvania State University
------------------------------------------------------------------------- */
#include <mpi.h>
#include <cstring>
#include <cstdlib>
#include "compute_pressure_bocs.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "modify.h"
#include "fix.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "kspace.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputePressureBocs::ComputePressureBocs(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
vptr(NULL), id_temp(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal compute pressure/bocs command");
if (igroup) error->all(FLERR,"Compute pressure/bocs must use group all");
scalar_flag = vector_flag = 1;
size_vector = 6;
extscalar = 0;
extvector = 0;
pressflag = 1;
timeflag = 1;
p_match_flag = 0;
phi_coeff = NULL;
// store temperature ID used by pressure computation
// insure it is valid for temperature computation
if (strcmp(arg[3],"NULL") == 0) id_temp = NULL;
else {
int n = strlen(arg[3]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[3]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0)
error->all(FLERR,"Could not find compute pressure/bocs temperature ID");
if (modify->compute[icompute]->tempflag == 0)
error->all(FLERR,"Compute pressure/bocs temperature ID does not "
"compute temperature");
}
// process optional args
if (narg == 4) {
keflag = 1;
pairflag = 1;
bondflag = angleflag = dihedralflag = improperflag = 1;
kspaceflag = fixflag = 1;
} else {
keflag = 0;
pairflag = 0;
bondflag = angleflag = dihedralflag = improperflag = 0;
kspaceflag = fixflag = 0;
int iarg = 4;
while (iarg < narg) {
if (strcmp(arg[iarg],"ke") == 0) keflag = 1;
else if (strcmp(arg[iarg],"pair") == 0) pairflag = 1;
else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1;
else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1;
else if (strcmp(arg[iarg],"dihedral") == 0) dihedralflag = 1;
else if (strcmp(arg[iarg],"improper") == 0) improperflag = 1;
else if (strcmp(arg[iarg],"kspace") == 0) kspaceflag = 1;
else if (strcmp(arg[iarg],"fix") == 0) fixflag = 1;
else if (strcmp(arg[iarg],"virial") == 0) {
pairflag = 1;
bondflag = angleflag = dihedralflag = improperflag = 1;
kspaceflag = fixflag = 1;
} else error->all(FLERR,"Illegal compute pressure/bocs command");
iarg++;
}
}
// error check
if (keflag && id_temp == NULL)
error->all(FLERR,"Compute pressure/bocs requires temperature ID "
"to include kinetic energy");
vector = new double[6];
nvirial = 0;
vptr = NULL;
}
/* ---------------------------------------------------------------------- */
ComputePressureBocs::~ComputePressureBocs()
{
delete [] id_temp;
delete [] vector;
delete [] vptr;
if (phi_coeff) free(phi_coeff);
}
/* ---------------------------------------------------------------------- */
void ComputePressureBocs::init()
{
boltz = force->boltz;
nktv2p = force->nktv2p;
dimension = domain->dimension;
// set temperature compute, must be done in init()
// fixes could have changed or compute_modify could have changed it
if (keflag) {
int icompute = modify->find_compute(id_temp);
if (icompute < 0)
error->all(FLERR,"Could not find compute pressure/bocs temperature ID");
temperature = modify->compute[icompute];
}
// detect contributions to virial
// vptr points to all virial[6] contributions
delete [] vptr;
nvirial = 0;
vptr = NULL;
if (pairflag && force->pair) nvirial++;
if (bondflag && atom->molecular && force->bond) nvirial++;
if (angleflag && atom->molecular && force->angle) nvirial++;
if (dihedralflag && atom->molecular && force->dihedral) nvirial++;
if (improperflag && atom->molecular && force->improper) nvirial++;
if (fixflag)
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->virial_flag) nvirial++;
if (nvirial) {
vptr = new double*[nvirial];
nvirial = 0;
if (pairflag && force->pair) vptr[nvirial++] = force->pair->virial;
if (bondflag && force->bond) vptr[nvirial++] = force->bond->virial;
if (angleflag && force->angle) vptr[nvirial++] = force->angle->virial;
if (dihedralflag && force->dihedral)
vptr[nvirial++] = force->dihedral->virial;
if (improperflag && force->improper)
vptr[nvirial++] = force->improper->virial;
if (fixflag)
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->virial_flag)
vptr[nvirial++] = modify->fix[i]->virial;
}
// flag Kspace contribution separately, since not summed across procs
if (kspaceflag && force->kspace) kspace_virial = force->kspace->virial;
else kspace_virial = NULL;
}
/* Extra functions added for BOCS */
/* ----------------------------------------------------------------------
Compute the pressure correction for the analytical basis set
------------------------------------------------------------------------- */
double ComputePressureBocs::get_cg_p_corr(int N_basis, double *phi_coeff,
int N_mol, double vavg, double vCG)
{
double correction = 0.0;
for (int i = 1; i <= N_basis; ++i)
{
correction -= phi_coeff[i-1] * ( N_mol * i / vavg ) *
pow( ( 1 / vavg ) * ( vCG - vavg ),i-1);
}
return correction;
}
/* ----------------------------------------------------------------------
Find the relevant index position if using a spline basis set
------------------------------------------------------------------------- */
double ComputePressureBocs::find_index(double * grid, double value)
{
int i;
double spacing = fabs(grid[1]-grid[0]);
int gridsize = spline_length;
for (i = 0; i < (gridsize-1); ++i)
{
if (value >= grid[i] && value <= grid[i+1]) { return i; }
}
if (value >= grid[i] && value <= (grid[i] + spacing)) { return i; }
for (int i = 0; i < gridsize; ++i)
{
fprintf(stderr, "grid %d: %f\n",i,grid[i]);
}
char * errmsg = (char *) calloc(100,sizeof(char));
sprintf(errmsg,"Value %f does not fall within spline grid.\n",value);
error->all(FLERR,errmsg);
exit(1);
}
/* ----------------------------------------------------------------------
Compute the pressure correction for a spline basis set
------------------------------------------------------------------------- */
double ComputePressureBocs::get_cg_p_corr(double ** grid, int basis_type,
double vCG)
{
int i = find_index(grid[0],vCG);
double correction, deltax = vCG - grid[0][i];
if (basis_type == 1)
{
correction = grid[1][i] + (deltax) *
( grid[1][i+1] - grid[1][i] ) / ( grid[0][i+1] - grid[0][i] );
}
else if (basis_type == 2)
{
correction = grid[1][i] + (grid[2][i] * deltax) +
(grid[3][i] * pow(deltax,2)) + (grid[4][i] * pow(deltax,3));
}
else
{
error->all(FLERR,"bad spline type passed to get_cg_p_corr()\n");
}
return correction;
}
/* ----------------------------------------------------------------------
send cg info from fix_bocs to compute_pressure_bocs for the analytical
basis set
------------------------------------------------------------------------- */
void ComputePressureBocs::send_cg_info(int basis_type, int sent_N_basis,
double *sent_phi_coeff, int sent_N_mol, double sent_vavg)
{
if (basis_type == 0) { p_basis_type = 0; }
else
{
error->all(FLERR,"Incorrect basis type passed to ComputePressureBocs\n");
}
p_match_flag = 1;
N_basis = sent_N_basis;
if (phi_coeff) free(phi_coeff);
phi_coeff = ((double *) calloc(N_basis, sizeof(double)) );
for (int i=0; i<N_basis; i++) { phi_coeff[i] = sent_phi_coeff[i]; }
N_mol = sent_N_mol;
vavg = sent_vavg;
}
/* ----------------------------------------------------------------------
send cg info from fix_bocs to compute_pressure_bocs for a spline basis
set
------------------------------------------------------------------------- */
void ComputePressureBocs::send_cg_info(int basis_type,
double ** in_splines, int gridsize)
{
if (basis_type == 1) { p_basis_type = 1; }
else if (basis_type == 2) { p_basis_type = 2; }
else
{
error->all(FLERR,"Incorrect basis type passed to ComputePressureBocs\n");
}
splines = in_splines;
spline_length = gridsize;
p_match_flag = 1;
}
/* End of new functions for BOCS */
/* ----------------------------------------------------------------------
compute total pressure, averaged over Pxx, Pyy, Pzz
------------------------------------------------------------------------- */
double ComputePressureBocs::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;
double volume, correction = 0;
if (keflag) {
if (temperature->invoked_scalar != update->ntimestep)
t = temperature->compute_scalar();
else t = temperature->scalar;
}
if (dimension == 3) {
inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
volume = (domain->xprd * domain->yprd * domain->zprd);
/* MRD NJD if block */
if ( p_basis_type == 0 )
{
correction = get_cg_p_corr(N_basis,phi_coeff,N_mol,vavg,volume);
}
else if ( p_basis_type == 1 || p_basis_type == 2 )
{
correction = get_cg_p_corr(splines, p_basis_type, volume);
}
virial_compute(3,3);
if (keflag)
scalar = (temperature->dof * boltz * t +
virial[0] + virial[1] + virial[2]) / 3.0 *
inv_volume * nktv2p + (correction);
else
scalar = (virial[0] + virial[1] + virial[2]) / 3.0 *
inv_volume * nktv2p + (correction);
} else {
if (p_match_flag)
{
error->all(FLERR,"Pressure matching not implemented in 2-d.\n");
exit(1);
} // The rest of this can probably be deleted.
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 ComputePressureBocs::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;
if (keflag) {
if (temperature->invoked_vector != update->ntimestep)
temperature->compute_vector();
ke_tensor = temperature->vector;
}
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;
}
}
}
/* ---------------------------------------------------------------------- */
void ComputePressureBocs::virial_compute(int n, int ndiag)
{
int i,j;
double v[6],*vcomponent;
for (i = 0; i < n; i++) v[i] = 0.0;
// sum contributions to virial from forces and fixes
for (j = 0; j < nvirial; j++) {
vcomponent = vptr[j];
for (i = 0; i < n; i++) v[i] += vcomponent[i];
}
// sum virial across procs
MPI_Allreduce(v,virial,n,MPI_DOUBLE,MPI_SUM,world);
// KSpace virial contribution is already summed across procs
if (kspace_virial)
for (i = 0; i < n; i++) virial[i] += kspace_virial[i];
// LJ long-range tail correction, only if pair contributions are included
if (force->pair && pairflag && force->pair->tail_flag)
for (i = 0; i < ndiag; i++) virial[i] += force->pair->ptail * inv_volume;
}
/* ---------------------------------------------------------------------- */
void ComputePressureBocs::reset_extra_compute_fix(const char *id_new)
{
delete [] id_temp;
int n = strlen(id_new) + 1;
id_temp = new char[n];
strcpy(id_temp,id_new);
}

View File

@ -0,0 +1,114 @@
/* -*- 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.
-------------------------------------------------------------------------
USER-BOCS written by: Nicholas J. H. Dunn and Michael R. DeLyser
from The Pennsylvania State University
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(PRESSURE/BOCS,ComputePressureBocs)
#else
#ifndef LMP_COMPUTE_PRESSURE_BOCS_H
#define LMP_COMPUTE_PRESSURE_BOCS_H
#include "compute.h"
namespace LAMMPS_NS {
// ComputePressure -> ComputePressureBocs MRD NJD
class ComputePressureBocs : public Compute {
public:
ComputePressureBocs(class LAMMPS *, int, char **);
virtual ~ComputePressureBocs();
virtual void init();
virtual double compute_scalar();
virtual void compute_vector();
void reset_extra_compute_fix(const char *);
double compute_cg_scalar();
double get_cg_p_corr(int, double *, int, double, double);
double get_cg_fluct(double, double);
void send_cg_info(int, int, double*, int, double);
void send_cg_info(int, double **, int);
double get_cg_p_corr(double **, int, double);
double find_index(double* , double);
protected:
double boltz,nktv2p,inv_volume;
int nvirial,dimension;
double **vptr;
double *kspace_virial;
Compute *temperature;
char *id_temp;
double virial[6];
int keflag,pairflag,bondflag,angleflag,dihedralflag,improperflag;
int fixflag,kspaceflag;
// NJD MRD
int p_basis_type;
int p_match_flag;
double vavg;
int N_mol;
int N_basis;
double *phi_coeff;
double ** splines;
int spline_length;
void virial_compute(int, int);
};
}
#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.
*/

2431
src/USER-BOCS/fix_bocs.cpp Normal file

File diff suppressed because it is too large Load Diff

302
src/USER-BOCS/fix_bocs.h Normal file
View File

@ -0,0 +1,302 @@
/* -*- 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.
-------------------------------------------------------------------------
USER-BOCS written by: Nicholas J. H. Dunn and Michael R. DeLyser
from The Pennsylvania State University
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(bocs,FixBocs)
#else
#ifndef LMP_FIX_BOCS_H
#define LMP_FIX_BOCS_H
#include "fix.h"
namespace LAMMPS_NS {
class FixBocs : public Fix {
public:
FixBocs(class LAMMPS *, int, char **); // MRD NJD
virtual ~FixBocs(); // MRD NJD
int setmask();
virtual void init();
virtual void setup(int);
virtual void initial_integrate(int);
virtual void final_integrate();
void initial_integrate_respa(int, int, int);
void final_integrate_respa(int, int);
virtual void pre_exchange();
double compute_scalar();
virtual double compute_vector(int);
void write_restart(FILE *);
virtual int pack_restart_data(double *); // pack restart data
virtual void restart(char *);
int modify_param(int, char **);
void reset_target(double);
void reset_dt();
virtual void *extract(const char*,int &);
double memory_usage();
protected:
int dimension,which;
double dtv,dtf,dthalf,dt4,dt8,dto;
double boltz,nktv2p,tdof;
double vol0; // reference volume
double t0; // reference temperature
// used for barostat mass
double t_start,t_stop;
double t_current,t_target,ke_target;
double t_freq;
int tstat_flag; // 1 if control T
int pstat_flag; // 1 if control P
int pstyle,pcouple,allremap;
int p_flag[6]; // 1 if control P on this dim, 0 if not
double p_start[6],p_stop[6];
double p_freq[6],p_target[6];
double omega[6],omega_dot[6];
double omega_mass[6];
double p_current[6];
double drag,tdrag_factor; // drag factor on particle thermostat
double pdrag_factor; // drag factor on barostat
int kspace_flag; // 1 if KSpace invoked, 0 if not
int nrigid; // number of rigid fixes
int dilate_group_bit; // mask for dilation group
int *rfix; // indices of rigid fixes
char *id_dilate; // group name to dilate
class Irregular *irregular; // for migrating atoms after box flips
// MRD NJD
int p_basis_type;
int p_match_flag;
double vavg;
int N_mol;
int N_p_match;
double *p_match_coeffs;
double ** splines;
int spline_length;
int nlevels_respa;
double *step_respa;
char *id_temp,*id_press;
class Compute *temperature,*pressure;
int tcomputeflag,pcomputeflag; // 1 = compute was created by fix
// 0 = created externally
double *eta,*eta_dot; // chain thermostat for particles
double *eta_dotdot;
double *eta_mass;
int mtchain; // length of chain
int mtchain_default_flag; // 1 = mtchain is default
double *etap; // chain thermostat for barostat
double *etap_dot;
double *etap_dotdot;
double *etap_mass;
int mpchain; // length of chain
int mtk_flag; // 0 if using Hoover barostat
int pdim; // number of barostatted dims
double p_freq_max; // maximum barostat frequency
double p_hydro; // hydrostatic target pressure
int nc_tchain,nc_pchain;
double factor_eta;
double sigma[6]; // scaled target stress
double fdev[6]; // deviatoric force on barostat
int deviatoric_flag; // 0 if target stress tensor is hydrostatic
double h0_inv[6]; // h_inv of reference (zero strain) box
int nreset_h0; // interval for resetting h0
double mtk_term1,mtk_term2; // Martyna-Tobias-Klein corrections
int eta_mass_flag; // 1 if eta_mass updated, 0 if not.
int omega_mass_flag; // 1 if omega_mass updated, 0 if not.
int etap_mass_flag; // 1 if etap_mass updated, 0 if not.
int dipole_flag; // 1 if dipole is updated, 0 if not.
int dlm_flag; // 1 if using the DLM rotational integrator, 0 if not
int scaleyz; // 1 if yz scaled with lz
int scalexz; // 1 if xz scaled with lz
int scalexy; // 1 if xy scaled with ly
int flipflag; // 1 if box flips are invoked as needed
int pre_exchange_flag; // set if pre_exchange needed for box flips
double fixedpoint[3]; // location of dilation fixed-point
void couple();
virtual void remap();
void nhc_temp_integrate();
void nhc_press_integrate();
int read_F_table(char *, int);
void build_cubic_splines(double **);
virtual void nve_x(); // may be overwritten by child classes
virtual void nve_v();
virtual void nh_v_press();
virtual void nh_v_temp();
virtual void compute_temp_target();
virtual int size_restart_global();
void compute_sigma();
void compute_deviatoric();
double compute_strain_energy();
void compute_press_target();
void nh_omega_dot();
};
}
#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: Target temperature for fix bocs cannot be 0.0
Self-explanatory.
E: Invalid fix bocs command for a 2d simulation
Cannot control z dimension in a 2d model.
E: Fix bocs dilate group ID does not exist
Self-explanatory.
E: Invalid fix bocs command pressure settings
If multiple dimensions are coupled, those dimensions must be
specified.
E: Cannot use fix bocs on a non-periodic dimension
When specifying a diagonal pressure component, the dimension must be
periodic.
E: Cannot use fix bocs on a 2nd non-periodic dimension
When specifying an off-diagonal pressure component, the 2nd of the two
dimensions must be periodic. E.g. if the xy component is specified,
then the y dimension must be periodic.
E: Cannot use fix bocs with yz scaling when z is non-periodic dimension
The 2nd dimension in the barostatted tilt factor must be periodic.
E: Cannot use fix bocs with xz scaling when z is non-periodic dimension
The 2nd dimension in the barostatted tilt factor must be periodic.
E: Cannot use fix bocs with xy scaling when y is non-periodic dimension
The 2nd dimension in the barostatted tilt factor must be periodic.
E: Cannot use fix bocs with both yz dynamics and yz scaling
Self-explanatory.
E: Cannot use fix bocs with both xz dynamics and xz scaling
Self-explanatory.
E: Cannot use fix bocs with both xy dynamics and xy scaling
Self-explanatory.
E: Can not specify Pxy/Pxz/Pyz in fix bocs with non-triclinic box
Only triclinic boxes can be used with off-diagonal pressure components.
See the region prism command for details.
E: Invalid fix bocs pressure settings
Settings for coupled dimensions must be the same.
E: Using update dipole flag requires atom style sphere
Self-explanatory.
E: Using update dipole flag requires atom attribute mu
Self-explanatory.
E: The dlm flag must be used with update dipole
Self-explanatory.
E: Fix bocs damping parameters must be > 0.0
Self-explanatory.
E: Cannot use fix npt and fix deform on same component of stress tensor
This would be changing the same box dimension twice.
E: Temperature ID for fix bocs does not exist
Self-explanatory.
E: Pressure ID for fix bocs does not exist
Self-explanatory.
E: Fix bocs has tilted box too far in one step - periodic cell is too far from equilibrium state
Self-explanatory. The change in the box tilt is too extreme
on a short timescale.
E: Could not find fix_modify temperature ID
The compute ID for computing temperature does not exist.
E: Fix_modify temperature ID does not compute temperature
The compute ID assigned to the fix must compute temperature.
W: Temperature for fix modify is not for group all
The temperature compute is being used with a pressure calculation
which does operate on group all, so this may be inconsistent.
E: Pressure ID for fix modify does not exist
Self-explanatory.
E: Could not find fix_modify pressure ID
The compute ID for computing pressure does not exist.
E: Fix_modify pressure ID does not compute pressure
The compute ID assigned to the fix must compute pressure.
*/