diff --git a/Documentation/source/guides/hma.rst b/Documentation/source/guides/hma.rst new file mode 100644 index 00000000..7c20e40d --- /dev/null +++ b/Documentation/source/guides/hma.rst @@ -0,0 +1,42 @@ +.. highlight:: none + +.. _ch:hma: + +Harmonically Mapped Averaging (HMA) +=================================== + +.. _sec:hma: + +Background +---------- + +Harmonically mapped averaging (HMA) is a technique that improves the precision +of measured themodynamic properites of crystalline systems. Beyond precision, +other benefits include smaller potential-truncation effects, finite-size +effects, faster equilibration and shorter decorrelation time. + +Limitations +----------- + +HMA can only be used in the NVT ensemble and when the atoms only vibrate around +their lattice sites (without diffusing). HMA can be used in Cassandra with any +pair potential, but not with Ewald summation. + +Truncation +---------- + +In order to supress property fluctuations due to atoms moving into and out of +the truncation radius, the potential in Cassandra is truncated based on the +lattice site distance instead of the distance between the atoms. This will lead +to a small difference in the measured properties, but the difference will +disappear in the limit of long truncation. + +Example +------- + +An simuation that uses HMA is included in `Examples/NVT/HMA/ `_ in the source tree. The example simulation runs a Lennard-Jones FCC crystal near its melting point. Energy and pressure are computed using both conventional and HMA and reported in the property output. + +References +---------- + +Sabry G. Moustafa, Andrew J. Schultz, and David A. Kofke, Very fast averaging of thermal properties of crystals by molecular simulation, `Phys. Rev. E [92], 043303 (2015) `_ diff --git a/Documentation/source/guides/input_files.rst b/Documentation/source/guides/input_files.rst index 15aa30cf..49b40db0 100644 --- a/Documentation/source/guides/input_files.rst +++ b/Documentation/source/guides/input_files.rst @@ -1353,6 +1353,8 @@ Additional supported keywords are: * ``nmols``: Number of molecules of each species * ``density``: Density of each species in Å\ :sup:`3` * ``mass_density``: Density of the system in kg/m\ :sup:`3` +* ``energy_HMA``: Energy of the system measured with :ref:`HMA ` +* ``pressure_HMA``: Pressure of the system measured with :ref:`HMA `, in bar. The harmonic pressure must be provided as a second argument on this line. For example, if you would like total energy, volume and pressure of a one box system to be written, you may specify the following: diff --git a/Documentation/source/index.rst b/Documentation/source/index.rst index 7e2cdade..7fd6f447 100644 --- a/Documentation/source/index.rst +++ b/Documentation/source/index.rst @@ -24,6 +24,7 @@ Cassandra Monte Carlo Software guides/input_files guides/utilities guides/solids + guides/hma .. toctree:: :hidden: diff --git a/Examples/NVT/HMA/README b/Examples/NVT/HMA/README new file mode 100644 index 00000000..8a31e420 --- /dev/null +++ b/Examples/NVT/HMA/README @@ -0,0 +1,19 @@ +This directory contains example HMA input files to run a 500-atom FCC +simulation using harmonically mapped averaging (HMA) with a Lennard-Jones +potential. + +More information about HMA is available in these publications: + +A. J. Schultz, D. A. Kofke, “Comprehensive high-precision high-accuracy +equation of state and coexistence properties for classical Lennard-Jones +crystals and low-temperature fluid phases”, J. Chem. Phys. 149, 204508 (2018) +https://doi.org/10.1063/1.5053714 + +S. G. Moustafa, A. J. Schultz, D. A. Kofke, “Harmonically Assisted Methods for +Computing the Free Energy of Classical Crystals by Molecular Simulation: A +Comparative Study”, J. Chem. Theory Comput. 13, 825-834 (2017) +https://doi.org/10.1021/acs.jctc.6b01082 + +S. G. Moustafa, A. J. Schultz, D. A. Kofke, “Very fast averaging of thermal +properties of crystals by molecular simulation”, Phys. Rev. E 92, 043303 (2015) +https://doi.org/10.1103/PhysRevE.92.043303 diff --git a/Examples/NVT/HMA/fcc500.xyz b/Examples/NVT/HMA/fcc500.xyz new file mode 100644 index 00000000..bd43c226 --- /dev/null +++ b/Examples/NVT/HMA/fcc500.xyz @@ -0,0 +1,502 @@ +500 +# lattice +A -3.571652366928449 -3.571652366928449 -3.571652366928449 +A -2.777951840944349 -2.777951840944349 -3.571652366928449 +A -2.777951840944349 -3.571652366928449 -2.777951840944349 +A -3.571652366928449 -2.777951840944349 -2.777951840944349 +A -3.571652366928449 -3.571652366928449 -1.984251314960249 +A -2.777951840944349 -2.777951840944349 -1.984251314960249 +A -2.777951840944349 -3.571652366928449 -1.190550788976150 +A -3.571652366928449 -2.777951840944349 -1.190550788976150 +A -3.571652366928449 -3.571652366928449 -0.396850262992050 +A -2.777951840944349 -2.777951840944349 -0.396850262992050 +A -2.777951840944349 -3.571652366928449 0.396850262992050 +A -3.571652366928449 -2.777951840944349 0.396850262992050 +A -3.571652366928449 -3.571652366928449 1.190550788976150 +A -2.777951840944349 -2.777951840944349 1.190550788976150 +A -2.777951840944349 -3.571652366928449 1.984251314960249 +A -3.571652366928449 -2.777951840944349 1.984251314960249 +A -3.571652366928449 -3.571652366928449 2.777951840944349 +A -2.777951840944349 -2.777951840944349 2.777951840944349 +A -2.777951840944349 -3.571652366928449 3.571652366928448 +A -3.571652366928449 -2.777951840944349 3.571652366928448 +A -3.571652366928449 -1.984251314960249 -3.571652366928449 +A -2.777951840944349 -1.190550788976150 -3.571652366928449 +A -2.777951840944349 -1.984251314960249 -2.777951840944349 +A -3.571652366928449 -1.190550788976150 -2.777951840944349 +A -3.571652366928449 -1.984251314960249 -1.984251314960249 +A -2.777951840944349 -1.190550788976150 -1.984251314960249 +A -2.777951840944349 -1.984251314960249 -1.190550788976150 +A -3.571652366928449 -1.190550788976150 -1.190550788976150 +A -3.571652366928449 -1.984251314960249 -0.396850262992050 +A -2.777951840944349 -1.190550788976150 -0.396850262992050 +A -2.777951840944349 -1.984251314960249 0.396850262992050 +A -3.571652366928449 -1.190550788976150 0.396850262992050 +A -3.571652366928449 -1.984251314960249 1.190550788976150 +A -2.777951840944349 -1.190550788976150 1.190550788976150 +A -2.777951840944349 -1.984251314960249 1.984251314960249 +A -3.571652366928449 -1.190550788976150 1.984251314960249 +A -3.571652366928449 -1.984251314960249 2.777951840944349 +A -2.777951840944349 -1.190550788976150 2.777951840944349 +A -2.777951840944349 -1.984251314960249 3.571652366928448 +A -3.571652366928449 -1.190550788976150 3.571652366928448 +A -3.571652366928449 -0.396850262992050 -3.571652366928449 +A -2.777951840944349 0.396850262992050 -3.571652366928449 +A -2.777951840944349 -0.396850262992050 -2.777951840944349 +A -3.571652366928449 0.396850262992050 -2.777951840944349 +A -3.571652366928449 -0.396850262992050 -1.984251314960249 +A -2.777951840944349 0.396850262992050 -1.984251314960249 +A -2.777951840944349 -0.396850262992050 -1.190550788976150 +A -3.571652366928449 0.396850262992050 -1.190550788976150 +A -3.571652366928449 -0.396850262992050 -0.396850262992050 +A -2.777951840944349 0.396850262992050 -0.396850262992050 +A -2.777951840944349 -0.396850262992050 0.396850262992050 +A -3.571652366928449 0.396850262992050 0.396850262992050 +A -3.571652366928449 -0.396850262992050 1.190550788976150 +A -2.777951840944349 0.396850262992050 1.190550788976150 +A -2.777951840944349 -0.396850262992050 1.984251314960249 +A -3.571652366928449 0.396850262992050 1.984251314960249 +A -3.571652366928449 -0.396850262992050 2.777951840944349 +A -2.777951840944349 0.396850262992050 2.777951840944349 +A -2.777951840944349 -0.396850262992050 3.571652366928448 +A -3.571652366928449 0.396850262992050 3.571652366928448 +A -3.571652366928449 1.190550788976150 -3.571652366928449 +A -2.777951840944349 1.984251314960249 -3.571652366928449 +A -2.777951840944349 1.190550788976150 -2.777951840944349 +A -3.571652366928449 1.984251314960249 -2.777951840944349 +A -3.571652366928449 1.190550788976150 -1.984251314960249 +A -2.777951840944349 1.984251314960249 -1.984251314960249 +A -2.777951840944349 1.190550788976150 -1.190550788976150 +A -3.571652366928449 1.984251314960249 -1.190550788976150 +A -3.571652366928449 1.190550788976150 -0.396850262992050 +A -2.777951840944349 1.984251314960249 -0.396850262992050 +A -2.777951840944349 1.190550788976150 0.396850262992050 +A -3.571652366928449 1.984251314960249 0.396850262992050 +A -3.571652366928449 1.190550788976150 1.190550788976150 +A -2.777951840944349 1.984251314960249 1.190550788976150 +A -2.777951840944349 1.190550788976150 1.984251314960249 +A -3.571652366928449 1.984251314960249 1.984251314960249 +A -3.571652366928449 1.190550788976150 2.777951840944349 +A -2.777951840944349 1.984251314960249 2.777951840944349 +A -2.777951840944349 1.190550788976150 3.571652366928448 +A -3.571652366928449 1.984251314960249 3.571652366928448 +A -3.571652366928449 2.777951840944349 -3.571652366928449 +A -2.777951840944349 3.571652366928448 -3.571652366928449 +A -2.777951840944349 2.777951840944349 -2.777951840944349 +A -3.571652366928449 3.571652366928448 -2.777951840944349 +A -3.571652366928449 2.777951840944349 -1.984251314960249 +A -2.777951840944349 3.571652366928448 -1.984251314960249 +A -2.777951840944349 2.777951840944349 -1.190550788976150 +A -3.571652366928449 3.571652366928448 -1.190550788976150 +A -3.571652366928449 2.777951840944349 -0.396850262992050 +A -2.777951840944349 3.571652366928448 -0.396850262992050 +A -2.777951840944349 2.777951840944349 0.396850262992050 +A -3.571652366928449 3.571652366928448 0.396850262992050 +A -3.571652366928449 2.777951840944349 1.190550788976150 +A -2.777951840944349 3.571652366928448 1.190550788976150 +A -2.777951840944349 2.777951840944349 1.984251314960249 +A -3.571652366928449 3.571652366928448 1.984251314960249 +A -3.571652366928449 2.777951840944349 2.777951840944349 +A -2.777951840944349 3.571652366928448 2.777951840944349 +A -2.777951840944349 2.777951840944349 3.571652366928448 +A -3.571652366928449 3.571652366928448 3.571652366928448 +A -1.984251314960249 -3.571652366928449 -3.571652366928449 +A -1.190550788976150 -2.777951840944349 -3.571652366928449 +A -1.190550788976150 -3.571652366928449 -2.777951840944349 +A -1.984251314960249 -2.777951840944349 -2.777951840944349 +A -1.984251314960249 -3.571652366928449 -1.984251314960249 +A -1.190550788976150 -2.777951840944349 -1.984251314960249 +A -1.190550788976150 -3.571652366928449 -1.190550788976150 +A -1.984251314960249 -2.777951840944349 -1.190550788976150 +A -1.984251314960249 -3.571652366928449 -0.396850262992050 +A -1.190550788976150 -2.777951840944349 -0.396850262992050 +A -1.190550788976150 -3.571652366928449 0.396850262992050 +A -1.984251314960249 -2.777951840944349 0.396850262992050 +A -1.984251314960249 -3.571652366928449 1.190550788976150 +A -1.190550788976150 -2.777951840944349 1.190550788976150 +A -1.190550788976150 -3.571652366928449 1.984251314960249 +A -1.984251314960249 -2.777951840944349 1.984251314960249 +A -1.984251314960249 -3.571652366928449 2.777951840944349 +A -1.190550788976150 -2.777951840944349 2.777951840944349 +A -1.190550788976150 -3.571652366928449 3.571652366928448 +A -1.984251314960249 -2.777951840944349 3.571652366928448 +A -1.984251314960249 -1.984251314960249 -3.571652366928449 +A -1.190550788976150 -1.190550788976150 -3.571652366928449 +A -1.190550788976150 -1.984251314960249 -2.777951840944349 +A -1.984251314960249 -1.190550788976150 -2.777951840944349 +A -1.984251314960249 -1.984251314960249 -1.984251314960249 +A -1.190550788976150 -1.190550788976150 -1.984251314960249 +A -1.190550788976150 -1.984251314960249 -1.190550788976150 +A -1.984251314960249 -1.190550788976150 -1.190550788976150 +A -1.984251314960249 -1.984251314960249 -0.396850262992050 +A -1.190550788976150 -1.190550788976150 -0.396850262992050 +A -1.190550788976150 -1.984251314960249 0.396850262992050 +A -1.984251314960249 -1.190550788976150 0.396850262992050 +A -1.984251314960249 -1.984251314960249 1.190550788976150 +A -1.190550788976150 -1.190550788976150 1.190550788976150 +A -1.190550788976150 -1.984251314960249 1.984251314960249 +A -1.984251314960249 -1.190550788976150 1.984251314960249 +A -1.984251314960249 -1.984251314960249 2.777951840944349 +A -1.190550788976150 -1.190550788976150 2.777951840944349 +A -1.190550788976150 -1.984251314960249 3.571652366928448 +A -1.984251314960249 -1.190550788976150 3.571652366928448 +A -1.984251314960249 -0.396850262992050 -3.571652366928449 +A -1.190550788976150 0.396850262992050 -3.571652366928449 +A -1.190550788976150 -0.396850262992050 -2.777951840944349 +A -1.984251314960249 0.396850262992050 -2.777951840944349 +A -1.984251314960249 -0.396850262992050 -1.984251314960249 +A -1.190550788976150 0.396850262992050 -1.984251314960249 +A -1.190550788976150 -0.396850262992050 -1.190550788976150 +A -1.984251314960249 0.396850262992050 -1.190550788976150 +A -1.984251314960249 -0.396850262992050 -0.396850262992050 +A -1.190550788976150 0.396850262992050 -0.396850262992050 +A -1.190550788976150 -0.396850262992050 0.396850262992050 +A -1.984251314960249 0.396850262992050 0.396850262992050 +A -1.984251314960249 -0.396850262992050 1.190550788976150 +A -1.190550788976150 0.396850262992050 1.190550788976150 +A -1.190550788976150 -0.396850262992050 1.984251314960249 +A -1.984251314960249 0.396850262992050 1.984251314960249 +A -1.984251314960249 -0.396850262992050 2.777951840944349 +A -1.190550788976150 0.396850262992050 2.777951840944349 +A -1.190550788976150 -0.396850262992050 3.571652366928448 +A -1.984251314960249 0.396850262992050 3.571652366928448 +A -1.984251314960249 1.190550788976150 -3.571652366928449 +A -1.190550788976150 1.984251314960249 -3.571652366928449 +A -1.190550788976150 1.190550788976150 -2.777951840944349 +A -1.984251314960249 1.984251314960249 -2.777951840944349 +A -1.984251314960249 1.190550788976150 -1.984251314960249 +A -1.190550788976150 1.984251314960249 -1.984251314960249 +A -1.190550788976150 1.190550788976150 -1.190550788976150 +A -1.984251314960249 1.984251314960249 -1.190550788976150 +A -1.984251314960249 1.190550788976150 -0.396850262992050 +A -1.190550788976150 1.984251314960249 -0.396850262992050 +A -1.190550788976150 1.190550788976150 0.396850262992050 +A -1.984251314960249 1.984251314960249 0.396850262992050 +A -1.984251314960249 1.190550788976150 1.190550788976150 +A -1.190550788976150 1.984251314960249 1.190550788976150 +A -1.190550788976150 1.190550788976150 1.984251314960249 +A -1.984251314960249 1.984251314960249 1.984251314960249 +A -1.984251314960249 1.190550788976150 2.777951840944349 +A -1.190550788976150 1.984251314960249 2.777951840944349 +A -1.190550788976150 1.190550788976150 3.571652366928448 +A -1.984251314960249 1.984251314960249 3.571652366928448 +A -1.984251314960249 2.777951840944349 -3.571652366928449 +A -1.190550788976150 3.571652366928448 -3.571652366928449 +A -1.190550788976150 2.777951840944349 -2.777951840944349 +A -1.984251314960249 3.571652366928448 -2.777951840944349 +A -1.984251314960249 2.777951840944349 -1.984251314960249 +A -1.190550788976150 3.571652366928448 -1.984251314960249 +A -1.190550788976150 2.777951840944349 -1.190550788976150 +A -1.984251314960249 3.571652366928448 -1.190550788976150 +A -1.984251314960249 2.777951840944349 -0.396850262992050 +A -1.190550788976150 3.571652366928448 -0.396850262992050 +A -1.190550788976150 2.777951840944349 0.396850262992050 +A -1.984251314960249 3.571652366928448 0.396850262992050 +A -1.984251314960249 2.777951840944349 1.190550788976150 +A -1.190550788976150 3.571652366928448 1.190550788976150 +A -1.190550788976150 2.777951840944349 1.984251314960249 +A -1.984251314960249 3.571652366928448 1.984251314960249 +A -1.984251314960249 2.777951840944349 2.777951840944349 +A -1.190550788976150 3.571652366928448 2.777951840944349 +A -1.190550788976150 2.777951840944349 3.571652366928448 +A -1.984251314960249 3.571652366928448 3.571652366928448 +A -0.396850262992050 -3.571652366928449 -3.571652366928449 +A 0.396850262992050 -2.777951840944349 -3.571652366928449 +A 0.396850262992050 -3.571652366928449 -2.777951840944349 +A -0.396850262992050 -2.777951840944349 -2.777951840944349 +A -0.396850262992050 -3.571652366928449 -1.984251314960249 +A 0.396850262992050 -2.777951840944349 -1.984251314960249 +A 0.396850262992050 -3.571652366928449 -1.190550788976150 +A -0.396850262992050 -2.777951840944349 -1.190550788976150 +A -0.396850262992050 -3.571652366928449 -0.396850262992050 +A 0.396850262992050 -2.777951840944349 -0.396850262992050 +A 0.396850262992050 -3.571652366928449 0.396850262992050 +A -0.396850262992050 -2.777951840944349 0.396850262992050 +A -0.396850262992050 -3.571652366928449 1.190550788976150 +A 0.396850262992050 -2.777951840944349 1.190550788976150 +A 0.396850262992050 -3.571652366928449 1.984251314960249 +A -0.396850262992050 -2.777951840944349 1.984251314960249 +A -0.396850262992050 -3.571652366928449 2.777951840944349 +A 0.396850262992050 -2.777951840944349 2.777951840944349 +A 0.396850262992050 -3.571652366928449 3.571652366928448 +A -0.396850262992050 -2.777951840944349 3.571652366928448 +A -0.396850262992050 -1.984251314960249 -3.571652366928449 +A 0.396850262992050 -1.190550788976150 -3.571652366928449 +A 0.396850262992050 -1.984251314960249 -2.777951840944349 +A -0.396850262992050 -1.190550788976150 -2.777951840944349 +A -0.396850262992050 -1.984251314960249 -1.984251314960249 +A 0.396850262992050 -1.190550788976150 -1.984251314960249 +A 0.396850262992050 -1.984251314960249 -1.190550788976150 +A -0.396850262992050 -1.190550788976150 -1.190550788976150 +A -0.396850262992050 -1.984251314960249 -0.396850262992050 +A 0.396850262992050 -1.190550788976150 -0.396850262992050 +A 0.396850262992050 -1.984251314960249 0.396850262992050 +A -0.396850262992050 -1.190550788976150 0.396850262992050 +A -0.396850262992050 -1.984251314960249 1.190550788976150 +A 0.396850262992050 -1.190550788976150 1.190550788976150 +A 0.396850262992050 -1.984251314960249 1.984251314960249 +A -0.396850262992050 -1.190550788976150 1.984251314960249 +A -0.396850262992050 -1.984251314960249 2.777951840944349 +A 0.396850262992050 -1.190550788976150 2.777951840944349 +A 0.396850262992050 -1.984251314960249 3.571652366928448 +A -0.396850262992050 -1.190550788976150 3.571652366928448 +A -0.396850262992050 -0.396850262992050 -3.571652366928449 +A 0.396850262992050 0.396850262992050 -3.571652366928449 +A 0.396850262992050 -0.396850262992050 -2.777951840944349 +A -0.396850262992050 0.396850262992050 -2.777951840944349 +A -0.396850262992050 -0.396850262992050 -1.984251314960249 +A 0.396850262992050 0.396850262992050 -1.984251314960249 +A 0.396850262992050 -0.396850262992050 -1.190550788976150 +A -0.396850262992050 0.396850262992050 -1.190550788976150 +A -0.396850262992050 -0.396850262992050 -0.396850262992050 +A 0.396850262992050 0.396850262992050 -0.396850262992050 +A 0.396850262992050 -0.396850262992050 0.396850262992050 +A -0.396850262992050 0.396850262992050 0.396850262992050 +A -0.396850262992050 -0.396850262992050 1.190550788976150 +A 0.396850262992050 0.396850262992050 1.190550788976150 +A 0.396850262992050 -0.396850262992050 1.984251314960249 +A -0.396850262992050 0.396850262992050 1.984251314960249 +A -0.396850262992050 -0.396850262992050 2.777951840944349 +A 0.396850262992050 0.396850262992050 2.777951840944349 +A 0.396850262992050 -0.396850262992050 3.571652366928448 +A -0.396850262992050 0.396850262992050 3.571652366928448 +A -0.396850262992050 1.190550788976150 -3.571652366928449 +A 0.396850262992050 1.984251314960249 -3.571652366928449 +A 0.396850262992050 1.190550788976150 -2.777951840944349 +A -0.396850262992050 1.984251314960249 -2.777951840944349 +A -0.396850262992050 1.190550788976150 -1.984251314960249 +A 0.396850262992050 1.984251314960249 -1.984251314960249 +A 0.396850262992050 1.190550788976150 -1.190550788976150 +A -0.396850262992050 1.984251314960249 -1.190550788976150 +A -0.396850262992050 1.190550788976150 -0.396850262992050 +A 0.396850262992050 1.984251314960249 -0.396850262992050 +A 0.396850262992050 1.190550788976150 0.396850262992050 +A -0.396850262992050 1.984251314960249 0.396850262992050 +A -0.396850262992050 1.190550788976150 1.190550788976150 +A 0.396850262992050 1.984251314960249 1.190550788976150 +A 0.396850262992050 1.190550788976150 1.984251314960249 +A -0.396850262992050 1.984251314960249 1.984251314960249 +A -0.396850262992050 1.190550788976150 2.777951840944349 +A 0.396850262992050 1.984251314960249 2.777951840944349 +A 0.396850262992050 1.190550788976150 3.571652366928448 +A -0.396850262992050 1.984251314960249 3.571652366928448 +A -0.396850262992050 2.777951840944349 -3.571652366928449 +A 0.396850262992050 3.571652366928448 -3.571652366928449 +A 0.396850262992050 2.777951840944349 -2.777951840944349 +A -0.396850262992050 3.571652366928448 -2.777951840944349 +A -0.396850262992050 2.777951840944349 -1.984251314960249 +A 0.396850262992050 3.571652366928448 -1.984251314960249 +A 0.396850262992050 2.777951840944349 -1.190550788976150 +A -0.396850262992050 3.571652366928448 -1.190550788976150 +A -0.396850262992050 2.777951840944349 -0.396850262992050 +A 0.396850262992050 3.571652366928448 -0.396850262992050 +A 0.396850262992050 2.777951840944349 0.396850262992050 +A -0.396850262992050 3.571652366928448 0.396850262992050 +A -0.396850262992050 2.777951840944349 1.190550788976150 +A 0.396850262992050 3.571652366928448 1.190550788976150 +A 0.396850262992050 2.777951840944349 1.984251314960249 +A -0.396850262992050 3.571652366928448 1.984251314960249 +A -0.396850262992050 2.777951840944349 2.777951840944349 +A 0.396850262992050 3.571652366928448 2.777951840944349 +A 0.396850262992050 2.777951840944349 3.571652366928448 +A -0.396850262992050 3.571652366928448 3.571652366928448 +A 1.190550788976150 -3.571652366928449 -3.571652366928449 +A 1.984251314960249 -2.777951840944349 -3.571652366928449 +A 1.984251314960249 -3.571652366928449 -2.777951840944349 +A 1.190550788976150 -2.777951840944349 -2.777951840944349 +A 1.190550788976150 -3.571652366928449 -1.984251314960249 +A 1.984251314960249 -2.777951840944349 -1.984251314960249 +A 1.984251314960249 -3.571652366928449 -1.190550788976150 +A 1.190550788976150 -2.777951840944349 -1.190550788976150 +A 1.190550788976150 -3.571652366928449 -0.396850262992050 +A 1.984251314960249 -2.777951840944349 -0.396850262992050 +A 1.984251314960249 -3.571652366928449 0.396850262992050 +A 1.190550788976150 -2.777951840944349 0.396850262992050 +A 1.190550788976150 -3.571652366928449 1.190550788976150 +A 1.984251314960249 -2.777951840944349 1.190550788976150 +A 1.984251314960249 -3.571652366928449 1.984251314960249 +A 1.190550788976150 -2.777951840944349 1.984251314960249 +A 1.190550788976150 -3.571652366928449 2.777951840944349 +A 1.984251314960249 -2.777951840944349 2.777951840944349 +A 1.984251314960249 -3.571652366928449 3.571652366928448 +A 1.190550788976150 -2.777951840944349 3.571652366928448 +A 1.190550788976150 -1.984251314960249 -3.571652366928449 +A 1.984251314960249 -1.190550788976150 -3.571652366928449 +A 1.984251314960249 -1.984251314960249 -2.777951840944349 +A 1.190550788976150 -1.190550788976150 -2.777951840944349 +A 1.190550788976150 -1.984251314960249 -1.984251314960249 +A 1.984251314960249 -1.190550788976150 -1.984251314960249 +A 1.984251314960249 -1.984251314960249 -1.190550788976150 +A 1.190550788976150 -1.190550788976150 -1.190550788976150 +A 1.190550788976150 -1.984251314960249 -0.396850262992050 +A 1.984251314960249 -1.190550788976150 -0.396850262992050 +A 1.984251314960249 -1.984251314960249 0.396850262992050 +A 1.190550788976150 -1.190550788976150 0.396850262992050 +A 1.190550788976150 -1.984251314960249 1.190550788976150 +A 1.984251314960249 -1.190550788976150 1.190550788976150 +A 1.984251314960249 -1.984251314960249 1.984251314960249 +A 1.190550788976150 -1.190550788976150 1.984251314960249 +A 1.190550788976150 -1.984251314960249 2.777951840944349 +A 1.984251314960249 -1.190550788976150 2.777951840944349 +A 1.984251314960249 -1.984251314960249 3.571652366928448 +A 1.190550788976150 -1.190550788976150 3.571652366928448 +A 1.190550788976150 -0.396850262992050 -3.571652366928449 +A 1.984251314960249 0.396850262992050 -3.571652366928449 +A 1.984251314960249 -0.396850262992050 -2.777951840944349 +A 1.190550788976150 0.396850262992050 -2.777951840944349 +A 1.190550788976150 -0.396850262992050 -1.984251314960249 +A 1.984251314960249 0.396850262992050 -1.984251314960249 +A 1.984251314960249 -0.396850262992050 -1.190550788976150 +A 1.190550788976150 0.396850262992050 -1.190550788976150 +A 1.190550788976150 -0.396850262992050 -0.396850262992050 +A 1.984251314960249 0.396850262992050 -0.396850262992050 +A 1.984251314960249 -0.396850262992050 0.396850262992050 +A 1.190550788976150 0.396850262992050 0.396850262992050 +A 1.190550788976150 -0.396850262992050 1.190550788976150 +A 1.984251314960249 0.396850262992050 1.190550788976150 +A 1.984251314960249 -0.396850262992050 1.984251314960249 +A 1.190550788976150 0.396850262992050 1.984251314960249 +A 1.190550788976150 -0.396850262992050 2.777951840944349 +A 1.984251314960249 0.396850262992050 2.777951840944349 +A 1.984251314960249 -0.396850262992050 3.571652366928448 +A 1.190550788976150 0.396850262992050 3.571652366928448 +A 1.190550788976150 1.190550788976150 -3.571652366928449 +A 1.984251314960249 1.984251314960249 -3.571652366928449 +A 1.984251314960249 1.190550788976150 -2.777951840944349 +A 1.190550788976150 1.984251314960249 -2.777951840944349 +A 1.190550788976150 1.190550788976150 -1.984251314960249 +A 1.984251314960249 1.984251314960249 -1.984251314960249 +A 1.984251314960249 1.190550788976150 -1.190550788976150 +A 1.190550788976150 1.984251314960249 -1.190550788976150 +A 1.190550788976150 1.190550788976150 -0.396850262992050 +A 1.984251314960249 1.984251314960249 -0.396850262992050 +A 1.984251314960249 1.190550788976150 0.396850262992050 +A 1.190550788976150 1.984251314960249 0.396850262992050 +A 1.190550788976150 1.190550788976150 1.190550788976150 +A 1.984251314960249 1.984251314960249 1.190550788976150 +A 1.984251314960249 1.190550788976150 1.984251314960249 +A 1.190550788976150 1.984251314960249 1.984251314960249 +A 1.190550788976150 1.190550788976150 2.777951840944349 +A 1.984251314960249 1.984251314960249 2.777951840944349 +A 1.984251314960249 1.190550788976150 3.571652366928448 +A 1.190550788976150 1.984251314960249 3.571652366928448 +A 1.190550788976150 2.777951840944349 -3.571652366928449 +A 1.984251314960249 3.571652366928448 -3.571652366928449 +A 1.984251314960249 2.777951840944349 -2.777951840944349 +A 1.190550788976150 3.571652366928448 -2.777951840944349 +A 1.190550788976150 2.777951840944349 -1.984251314960249 +A 1.984251314960249 3.571652366928448 -1.984251314960249 +A 1.984251314960249 2.777951840944349 -1.190550788976150 +A 1.190550788976150 3.571652366928448 -1.190550788976150 +A 1.190550788976150 2.777951840944349 -0.396850262992050 +A 1.984251314960249 3.571652366928448 -0.396850262992050 +A 1.984251314960249 2.777951840944349 0.396850262992050 +A 1.190550788976150 3.571652366928448 0.396850262992050 +A 1.190550788976150 2.777951840944349 1.190550788976150 +A 1.984251314960249 3.571652366928448 1.190550788976150 +A 1.984251314960249 2.777951840944349 1.984251314960249 +A 1.190550788976150 3.571652366928448 1.984251314960249 +A 1.190550788976150 2.777951840944349 2.777951840944349 +A 1.984251314960249 3.571652366928448 2.777951840944349 +A 1.984251314960249 2.777951840944349 3.571652366928448 +A 1.190550788976150 3.571652366928448 3.571652366928448 +A 2.777951840944349 -3.571652366928449 -3.571652366928449 +A 3.571652366928448 -2.777951840944349 -3.571652366928449 +A 3.571652366928448 -3.571652366928449 -2.777951840944349 +A 2.777951840944349 -2.777951840944349 -2.777951840944349 +A 2.777951840944349 -3.571652366928449 -1.984251314960249 +A 3.571652366928448 -2.777951840944349 -1.984251314960249 +A 3.571652366928448 -3.571652366928449 -1.190550788976150 +A 2.777951840944349 -2.777951840944349 -1.190550788976150 +A 2.777951840944349 -3.571652366928449 -0.396850262992050 +A 3.571652366928448 -2.777951840944349 -0.396850262992050 +A 3.571652366928448 -3.571652366928449 0.396850262992050 +A 2.777951840944349 -2.777951840944349 0.396850262992050 +A 2.777951840944349 -3.571652366928449 1.190550788976150 +A 3.571652366928448 -2.777951840944349 1.190550788976150 +A 3.571652366928448 -3.571652366928449 1.984251314960249 +A 2.777951840944349 -2.777951840944349 1.984251314960249 +A 2.777951840944349 -3.571652366928449 2.777951840944349 +A 3.571652366928448 -2.777951840944349 2.777951840944349 +A 3.571652366928448 -3.571652366928449 3.571652366928448 +A 2.777951840944349 -2.777951840944349 3.571652366928448 +A 2.777951840944349 -1.984251314960249 -3.571652366928449 +A 3.571652366928448 -1.190550788976150 -3.571652366928449 +A 3.571652366928448 -1.984251314960249 -2.777951840944349 +A 2.777951840944349 -1.190550788976150 -2.777951840944349 +A 2.777951840944349 -1.984251314960249 -1.984251314960249 +A 3.571652366928448 -1.190550788976150 -1.984251314960249 +A 3.571652366928448 -1.984251314960249 -1.190550788976150 +A 2.777951840944349 -1.190550788976150 -1.190550788976150 +A 2.777951840944349 -1.984251314960249 -0.396850262992050 +A 3.571652366928448 -1.190550788976150 -0.396850262992050 +A 3.571652366928448 -1.984251314960249 0.396850262992050 +A 2.777951840944349 -1.190550788976150 0.396850262992050 +A 2.777951840944349 -1.984251314960249 1.190550788976150 +A 3.571652366928448 -1.190550788976150 1.190550788976150 +A 3.571652366928448 -1.984251314960249 1.984251314960249 +A 2.777951840944349 -1.190550788976150 1.984251314960249 +A 2.777951840944349 -1.984251314960249 2.777951840944349 +A 3.571652366928448 -1.190550788976150 2.777951840944349 +A 3.571652366928448 -1.984251314960249 3.571652366928448 +A 2.777951840944349 -1.190550788976150 3.571652366928448 +A 2.777951840944349 -0.396850262992050 -3.571652366928449 +A 3.571652366928448 0.396850262992050 -3.571652366928449 +A 3.571652366928448 -0.396850262992050 -2.777951840944349 +A 2.777951840944349 0.396850262992050 -2.777951840944349 +A 2.777951840944349 -0.396850262992050 -1.984251314960249 +A 3.571652366928448 0.396850262992050 -1.984251314960249 +A 3.571652366928448 -0.396850262992050 -1.190550788976150 +A 2.777951840944349 0.396850262992050 -1.190550788976150 +A 2.777951840944349 -0.396850262992050 -0.396850262992050 +A 3.571652366928448 0.396850262992050 -0.396850262992050 +A 3.571652366928448 -0.396850262992050 0.396850262992050 +A 2.777951840944349 0.396850262992050 0.396850262992050 +A 2.777951840944349 -0.396850262992050 1.190550788976150 +A 3.571652366928448 0.396850262992050 1.190550788976150 +A 3.571652366928448 -0.396850262992050 1.984251314960249 +A 2.777951840944349 0.396850262992050 1.984251314960249 +A 2.777951840944349 -0.396850262992050 2.777951840944349 +A 3.571652366928448 0.396850262992050 2.777951840944349 +A 3.571652366928448 -0.396850262992050 3.571652366928448 +A 2.777951840944349 0.396850262992050 3.571652366928448 +A 2.777951840944349 1.190550788976150 -3.571652366928449 +A 3.571652366928448 1.984251314960249 -3.571652366928449 +A 3.571652366928448 1.190550788976150 -2.777951840944349 +A 2.777951840944349 1.984251314960249 -2.777951840944349 +A 2.777951840944349 1.190550788976150 -1.984251314960249 +A 3.571652366928448 1.984251314960249 -1.984251314960249 +A 3.571652366928448 1.190550788976150 -1.190550788976150 +A 2.777951840944349 1.984251314960249 -1.190550788976150 +A 2.777951840944349 1.190550788976150 -0.396850262992050 +A 3.571652366928448 1.984251314960249 -0.396850262992050 +A 3.571652366928448 1.190550788976150 0.396850262992050 +A 2.777951840944349 1.984251314960249 0.396850262992050 +A 2.777951840944349 1.190550788976150 1.190550788976150 +A 3.571652366928448 1.984251314960249 1.190550788976150 +A 3.571652366928448 1.190550788976150 1.984251314960249 +A 2.777951840944349 1.984251314960249 1.984251314960249 +A 2.777951840944349 1.190550788976150 2.777951840944349 +A 3.571652366928448 1.984251314960249 2.777951840944349 +A 3.571652366928448 1.190550788976150 3.571652366928448 +A 2.777951840944349 1.984251314960249 3.571652366928448 +A 2.777951840944349 2.777951840944349 -3.571652366928449 +A 3.571652366928448 3.571652366928448 -3.571652366928449 +A 3.571652366928448 2.777951840944349 -2.777951840944349 +A 2.777951840944349 3.571652366928448 -2.777951840944349 +A 2.777951840944349 2.777951840944349 -1.984251314960249 +A 3.571652366928448 3.571652366928448 -1.984251314960249 +A 3.571652366928448 2.777951840944349 -1.190550788976150 +A 2.777951840944349 3.571652366928448 -1.190550788976150 +A 2.777951840944349 2.777951840944349 -0.396850262992050 +A 3.571652366928448 3.571652366928448 -0.396850262992050 +A 3.571652366928448 2.777951840944349 0.396850262992050 +A 2.777951840944349 3.571652366928448 0.396850262992050 +A 2.777951840944349 2.777951840944349 1.190550788976150 +A 3.571652366928448 3.571652366928448 1.190550788976150 +A 3.571652366928448 2.777951840944349 1.984251314960249 +A 2.777951840944349 3.571652366928448 1.984251314960249 +A 2.777951840944349 2.777951840944349 2.777951840944349 +A 3.571652366928448 3.571652366928448 2.777951840944349 +A 3.571652366928448 2.777951840944349 3.571652366928448 +A 2.777951840944349 3.571652366928448 3.571652366928448 diff --git a/Examples/NVT/HMA/inp500_1.0 b/Examples/NVT/HMA/inp500_1.0 new file mode 100644 index 00000000..60396c08 --- /dev/null +++ b/Examples/NVT/HMA/inp500_1.0 @@ -0,0 +1,79 @@ +# Run_Name +out500 +!------------------------------------------------------------------------------ + +# Sim_Type +nvt +!------------------------------------------------------------------------------ + +# Nbr_Species +1 +!------------------------------------------------------------------------------ + +# Molecule_Files +lj.mcf 500 +!------------------------------------------------------------------------------ + +# Box_Info +1 +cubic +7.937005259840997 +!------------------------------------------------------------------------------ + +# VDW_Style +lj cut 3. +!---------------------------------------------------------one line for each box + +# Temperature_Info +1. +!---------------------------------------------------one entry/line for each box + +# Move_Probability_Info + +# Prob_Translation +100. +0.1 + +# Done_Probability_Info +!---------------------one line for each box and one entry/line for each species + +# Start_Type +read_config 500 fcc500.xyz +!---------------------one line for each species and one entry/line for each box + +# Run_Type +production 50000 100 +!------------------------------------------------------------------------------ + +# Simulation_Length_Info +units sweeps +prop_freq 10 +coord_freq 200 +run 2000 +!------------------------------------------------------------------------------ + +# Property_Info 1 +energy_total +energy_HMA +pressure +pressure_HMA 1311.62 +!------------------------------------------------------------------------------ + +# CBMC_Info +kappa_ins 12 +rcut_cbmc 6.5 +!------------------------------------------------------------------------------ + +# Rcutoff_Low +0.8 +!------------------------------------------------------------------------------ + +# Pair_Energy +true +!------------------------------------------------------------------------------ + +# Seed_Info +9799807 23445206 +!------------------------------------------------------------------------------ + +END diff --git a/Examples/NVT/HMA/lj.mcf b/Examples/NVT/HMA/lj.mcf new file mode 100644 index 00000000..b004374f --- /dev/null +++ b/Examples/NVT/HMA/lj.mcf @@ -0,0 +1,34 @@ +!******************************************************************************* +!Molecular connectivity file for n2.pdb +!******************************************************************************* +!Atom Format +!index type element mass charge vdw_type parameters +!vdw_type="LJ", parms=epsilon sigma +!vdw_type="Mie", parms=epsilon sigma repulsion_exponent dispersion_exponent + +# Atom_Info +1 +1 A A 1.000 0.000 LJ 1.000 1.000 + +!Bond Format +!index i j type parameters +!type="fixed", parms=bondLength + +# Bond_Info +0 + +!Angle Format +!index i j k type parameters +!type="fixed", parms=equilibrium_angle +!type="harmonic", parms=force_constant equilibrium_angle + +# Angle_Info +0 + +# Dihedral_Info +0 + +# Improper_Info +0 + +END diff --git a/Src/accumulate.f90 b/Src/accumulate.f90 index 6b33926c..8413ed32 100755 --- a/Src/accumulate.f90 +++ b/Src/accumulate.f90 @@ -115,6 +115,16 @@ SUBROUTINE Accumulate(ibox) mass_density = mass_density / box_list(ibox)%volume ac_mass_density(ibox,iblock) = ac_mass_density(ibox,iblock) + mass_density + IF (need_HMA) THEN + IF (energy_HMA(ibox)%last_calc /= i_mcstep) THEN + energy_HMA(ibox)%last_calc = i_mcstep + pressure_HMA(ibox)%last_calc = i_mcstep + CALL Compute_HMA(ibox) + END IF + ac_energy_HMA(ibox,iblock) = ac_energy_HMA(ibox,iblock) + energy_HMA(ibox)%total + ac_pressure_HMA(ibox,iblock) = ac_pressure_HMA(ibox,iblock) + pressure_HMA(ibox)%total + END IF + IF (MOD(i_mcstep,block_avg_freq) == 0) THEN ac_energy(ibox,iblock)%total = ac_energy(ibox,iblock)%total / data_points_per_block ac_energy(ibox,iblock)%intra = ac_energy(ibox,iblock)%intra / data_points_per_block @@ -138,5 +148,6 @@ SUBROUTINE Accumulate(ibox) ac_density(is,ibox,iblock) = ac_density(is,ibox,iblock) / data_points_per_block ac_nmols(is,ibox,iblock) = ac_nmols(is,ibox,iblock) / data_points_per_block END DO + ac_energy_HMA(ibox,iblock) = ac_energy_HMA(ibox,iblock) / data_points_per_block END IF END SUBROUTINE Accumulate diff --git a/Src/energy_routines.f90 b/Src/energy_routines.f90 index 2f94e7cf..8c0b378a 100755 --- a/Src/energy_routines.f90 +++ b/Src/energy_routines.f90 @@ -134,6 +134,11 @@ MODULE Energy_Routines ! ! Check_AtomPair_Cutoff: ! + ! Check_AtomPair_Lattice_Cutoff: Uses lattice sites to determine truncation + ! instead of molecules (or atoms) so that atoms do not + ! fluctuate in and out of the truncation during the + ! simulation. + ! ! Get_Molecule_Energy: Computes the intra- and inter-molecular energy of ! a given molecule interacting with all other molecules. ! @@ -654,7 +659,13 @@ SUBROUTINE Compute_Atom_Nonbond_Energy(this_atom,this_molecule,this_species, & ENDIF - CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box) + IF (need_HMA) THEN + CALL Check_AtomPair_Lattice_Cutoff(atom_list(ia,im,is), & + atom_list(this_atom,this_molecule,this_species), & + get_vdw,get_qq,this_box) + ELSE + CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box) + END IF ! Compute vdw and q-q energy using if required IF (get_vdw .OR. get_qq) THEN @@ -1121,7 +1132,12 @@ SUBROUTINE Compute_MoleculePair_Energy(im,is,jm,js,this_box, & END IF ! Now figure out what needs to be computed, then call pair_energy - CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box) + IF (need_HMA) THEN + CALL Check_AtomPair_Lattice_Cutoff(atom_list(ia,im,is), & + atom_list(ja,jm,js),get_vdw,get_qq,this_box) + ELSE + CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box) + END IF IF(cbmc_flag .AND. (.NOT. species_list(is)%L_Coul_CBMC)) THEN get_qq=.FALSE. @@ -2301,6 +2317,11 @@ SUBROUTINE Check_MoleculePair_Cutoff(im_1,is_1,im_2,is_2,get_interaction, & get_interaction = .TRUE. + IF (need_HMA) THEN + ! HMA uses an atomic cutoff on lattice sites + RETURN + END IF + ! Figure out the box to be used later. this_box = molecule_list(im_1,is_1)%which_box @@ -2446,6 +2467,36 @@ SUBROUTINE Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box) END SUBROUTINE Check_AtomPair_Cutoff +!******************************************************************************* + SUBROUTINE Check_AtomPair_Lattice_Cutoff(atom1,atom2,get_vdw,get_qq,this_box) + !*************************************************************************** + ! The subroutine determines if the interaction between atom1 and atom2 + ! should be truncated or not based on the distance between their lattice + ! sites. This routine computes the lattice distance and uses + ! Check_AtomPair_Cutoff to decide on the actual truncation. + !*************************************************************************** + TYPE(Atom_Class) :: atom1, atom2 + INTEGER :: this_box + LOGICAL :: get_vdw, get_qq + + REAL(DP) :: rxijp, ryijp, rzijp, rxij, ryij, rzij, rijsq + + rxijp = atom1%rxp_init - atom2%rxp_init + ryijp = atom1%ryp_init - atom2%ryp_init + rzijp = atom1%rzp_init - atom2%rzp_init + + ! Now get the minimum image separation + CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp,rxij,ryij,rzij) + + rijsq = rxij*rxij + ryij*ryij + rzij*rzij + + ! Now figure out what needs to be computed, then call pair_energy + CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box) + + RETURN + + END SUBROUTINE Check_AtomPair_Lattice_Cutoff + SUBROUTINE Compute_System_Total_Force(this_box) !**************************************************************************** @@ -2481,6 +2532,9 @@ SUBROUTINE Compute_System_Total_Force(this_box) W_tensor_charge(:,:,this_box) = 0.0_DP W_tensor_recip(:,:,this_box) = 0.0_DP W_tensor_elec(:,:,this_box) = 0.0_DP + IF (need_HMA) THEN + energy_HMA(this_box)%sum_Fdr = 0.0_DP + END IF DO is = 1, nspecies imLOOP1: DO im_1 = 1, nmols(is,this_box) @@ -2623,11 +2677,24 @@ SUBROUTINE Compute_MoleculePair_Force(im,is,jm,js,this_box,tens_vdw,tens_charge, rijsq = rxij*rxij + ryij*ryij + rzij*rzij ! Now figure out what needs to be computed, then call pair_energy - CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box) + IF (need_HMA) THEN + CALL Check_AtomPair_Lattice_Cutoff(atom_list(ia,im,is), & + atom_list(ja,jm,js),get_vdw,get_qq,this_box) + ELSE + CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box) + END IF ! Compute vdw and q-q energy using if required IF (get_vdw .OR. get_qq) THEN + IF (need_HMA) THEN + ! molecule-molecule vector was not computed before, but we + ! know the molecules have only one atom each + rabx = rxij + raby = ryij + rabz = rzij + END IF + CALL Compute_AtomPair_Force(rijsq,is,im,ia,js,jm,ja,& get_vdw,get_qq,Wij_vdw,Wij_qq) @@ -2685,7 +2752,7 @@ SUBROUTINE Compute_AtomPair_Force & ! Called by: Compute_System_Total_Force !----------------------------------------------------------------------------- ! Passed to - REAL(DP) :: rxij,ryij,rzij,rijsq + REAL(DP) :: rijsq INTEGER :: is,im,ia,js,jm,ja LOGICAL :: get_vdw,get_qq @@ -2706,6 +2773,11 @@ SUBROUTINE Compute_AtomPair_Force & ! Coulomb potential REAL(DP) :: qi, qj, erfc_val, prefactor REAL(DP) :: ewald_constant, exp_const + ! HMA + REAL(DP) :: dr_dot_product + REAL(DP) :: drixp, driyp, drizp, drix, driy, driz + REAL(DP) :: drjxp, drjyp, drjzp, drjx, drjy, drjz + REAL(DP) :: rxijp, ryijp, rzijp, rxij, ryij, rzij Wij_vdw = 0.0_DP Wij_qq = 0.0_DP @@ -2817,6 +2889,34 @@ SUBROUTINE Compute_AtomPair_Force & ENDIF qq_calc ENDIF ExistCheck + + HMA_calc: IF (need_HMA) THEN + drixp = atom_list(ia,im,is)%rxp - atom_list(ia,im,is)%rxp_init + driyp = atom_list(ia,im,is)%ryp - atom_list(ia,im,is)%ryp_init + drizp = atom_list(ia,im,is)%rzp - atom_list(ia,im,is)%rzp_init + CALL Minimum_Image_Separation(ibox,drixp,driyp,drizp, & + drix,driy,driz) + + drjxp = atom_list(ja,jm,js)%rxp - atom_list(ja,jm,js)%rxp_init + drjyp = atom_list(ja,jm,js)%ryp - atom_list(ja,jm,js)%ryp_init + drjzp = atom_list(ja,jm,js)%rzp - atom_list(ja,jm,js)%rzp_init + CALL Minimum_Image_Separation(ibox,drjxp,drjyp,drjzp, & + drjx,drjy,drjz) + + rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp + ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp + rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp + + ! Now get the minimum image separation + CALL Minimum_Image_Separation(ibox,rxijp,ryijp,rzijp, & + rxij,ryij,rzij) + + dr_dot_product = rxij*(drix-drjx) & + + ryij*(driy-drjy) & + + rzij*(driz-drjz) + energy_HMA(ibox)%sum_Fdr = energy_HMA(ibox)%sum_Fdr & + + (Wij_qq+Wij_vdw) * dr_dot_product / rijsq + ENDIF HMA_calc !------------------------------------------------------------------------------ CONTAINS diff --git a/Src/global_variables.f90 b/Src/global_variables.f90 index 8d8322bb..52e89c88 100755 --- a/Src/global_variables.f90 +++ b/Src/global_variables.f90 @@ -271,6 +271,13 @@ MODULE Global_Variables TYPE(Pressure_Class), DIMENSION(:), ALLOCATABLE, TARGET :: pressure LOGICAL :: need_pressure + !********************************************************************************** + ! HMA properties + + TYPE(Energy_HMA_Class), DIMENSION(:), ALLOCATABLE, TARGET :: energy_HMA + TYPE(Pressure_HMA_Class), DIMENSION(:), ALLOCATABLE, TARGET :: pressure_HMA + LOGICAL :: need_HMA + ! ********************************************************************************** ! system size integers used in memory allocation. ! Number of species, molecules, atoms, bonds, angles, dihedrals and impropers should @@ -423,6 +430,7 @@ MODULE Global_Variables ! will have dimensions of nbr_boxes REAL(DP), DIMENSION(:,:),ALLOCATABLE,TARGET :: ac_volume, ac_enthalpy, ac_pressure, ac_mass_density + REAL(DP), DIMENSION(:,:),ALLOCATABLE,TARGET :: ac_energy_HMA, ac_pressure_HMA ! will have dimension of (nspecies,nbr_boxes) REAL(DP), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ac_density, ac_nmols diff --git a/Src/initialize.f90 b/Src/initialize.f90 index cd3712fa..73ad91e1 100755 --- a/Src/initialize.f90 +++ b/Src/initialize.f90 @@ -135,6 +135,7 @@ SUBROUTINE Init_Accumulators ALLOCATE(ac_pressure(nbr_boxes,nbr_blocks)) ALLOCATE(ac_enthalpy(nbr_boxes,nbr_blocks)) ALLOCATE(ac_mass_density(nbr_boxes,nbr_blocks)) + ALLOCATE(ac_energy_HMA(nbr_boxes,nbr_blocks)) ac_energy(:,:)%total = 0.0_DP ac_energy(:,:)%intra = 0.0_DP @@ -163,6 +164,7 @@ SUBROUTINE Init_Accumulators ac_pressure = 0.0_DP ac_enthalpy = 0.0_DP ac_mass_density = 0.0_DP + ac_energy_HMA = 0.0_DP END IF END SUBROUTINE Init_Accumulators diff --git a/Src/input_routines.f90 b/Src/input_routines.f90 index a3494c9d..c200290a 100755 --- a/Src/input_routines.f90 +++ b/Src/input_routines.f90 @@ -5983,6 +5983,39 @@ SUBROUTINE Get_Property_Info ELSE IF (line_array(1) == 'enthalpy' .OR. line_array(1) == 'Enthalpy') THEN nbr_properties = nbr_properties + 1 IF (int_sim_type /= sim_npt .AND. int_sim_type /= sim_gemc_npt) need_pressure = .TRUE. + + ELSE IF (line_array(1) == 'energy_HMA' .OR. line_array(1) == 'Energy_HMA' .OR. & + line_array(1) == 'pressure_HMA' .OR. line_array(1) == 'Pressure_HMA') THEN + + IF (.NOT. need_HMA) THEN + ALLOCATE(energy_HMA(nbr_boxes)) + ALLOCATE(pressure_HMA(nbr_boxes)) + + energy_HMA(:)%lattice = 0.0_DP + energy_HMA(:)%harmonic = 0.0_DP + energy_HMA(:)%anharmonic = 0.0_DP + energy_HMA(:)%total = 0.0_DP + + pressure_HMA(:)%lattice = 0.0_DP + pressure_HMA(:)%harmonic = 0.0_DP + pressure_HMA(:)%anharmonic = 0.0_DP + pressure_HMA(:)%total = 0.0_DP + END IF + + need_HMA = .TRUE. + + IF (line_array(1) == 'pressure_HMA' .OR. line_array(1) == 'pressure_HMA') THEN + need_pressure = .TRUE. + ! harmonic pressure is the second value on the line + pressure_HMA(this_box)%harmonic = String_To_Double(line_array(2)) + WRITE(logunit,'(A,X,I1,X,A,X,F9.3,X,A)',ADVANCE='NO') 'HMA harmonic pressure of box', this_box, 'is', pressure_HMA(this_box)%harmonic, 'bar' + + ! convert pressure into atomic units + pressure_HMA(this_box)%harmonic = pressure_HMA(this_box)%harmonic / atomic_to_bar + WRITE(logunit,'(X,A,X,E13.6,X,A)') '=', pressure_HMA(this_box)%harmonic, 'amu / (A ps^2)' + END IF + + nbr_properties = nbr_properties + 1 ELSE ! this is a property for the system nbr_properties = nbr_properties + 1 @@ -6161,6 +6194,12 @@ SUBROUTINE Get_Property_Info ELSE IF (line_array(1) == 'mass_density' .OR. line_array(1) == 'Mass_Density') THEN nbr_properties = nbr_properties + 1 prop_output(nbr_properties,nbr_prop_files(this_box),this_box) = 'Mass_Density' + ELSE IF (line_array(1) == 'energy_HMA' .OR. line_array(1) == 'Energy_HMA') THEN + nbr_properties = nbr_properties + 1 + prop_output(nbr_properties,nbr_prop_files(this_box),this_box) = 'Energy_HMA' + ELSE IF (line_array(1) == 'pressure_HMA' .OR. line_array(1) == 'Pressure_HMA') THEN + nbr_properties = nbr_properties + 1 + prop_output(nbr_properties,nbr_prop_files(this_box),this_box) = 'Pressure_HMA' ELSE err_msg = '' err_msg(1) = 'Keyword "' // TRIM(line_array(1)) // '" on line ' // & @@ -6170,7 +6209,7 @@ SUBROUTINE Get_Property_Info err_msg(4) = ' energy_angle, energy_dihedral, energy_improper, energy_intravdw, energy_intraq' err_msg(5) = ' energy_intervdw, energy_interq, energy_lrc, energy_recip, energy_self,' err_msg(6) = ' enthalpy, pressure, pressure_xx, pressure_yy, pressure_zz, volume, density,' - err_msg(7) = ' nmols, mass_density' + err_msg(7) = ' nmols, mass_density, energy_HMA, pressure_HMA' CALL Clean_Abort(err_msg,'Get_Property_Info') END IF @@ -6243,7 +6282,7 @@ SUBROUTINE Get_Property_Info END IF - IF (need_pressure) THEN + IF (need_pressure .OR. need_HMA) THEN ALLOCATE(W_tensor_charge(3,3,nbr_boxes) , W_tensor_recip(3,3,nbr_boxes)) ALLOCATE(W_tensor_vdw(3,3,nbr_boxes) , W_tensor_total(3,3,nbr_boxes)) ALLOCATE(W_tensor_elec(3,3,nbr_boxes), pressure_tensor(3,3,nbr_boxes)) @@ -6702,5 +6741,4 @@ SUBROUTINE Get_Verbosity_Info WRITE(logunit,'(A80)') '********************************************************************************' END SUBROUTINE Get_Verbosity_Info - END MODULE Input_Routines diff --git a/Src/main.f90 b/Src/main.f90 index 1e5f38f4..d40fada5 100755 --- a/Src/main.f90 +++ b/Src/main.f90 @@ -378,6 +378,11 @@ PROGRAM Main WRITE(logunit,'(A80)') '********************************************************************************' DO ibox = 1,nbr_boxes CALL Check_System_Energy(ibox,.FALSE.) + IF (need_HMA) THEN + CALL Compute_Pressure(ibox) + energy_HMA(ibox)%lattice = energy(ibox)%total + pressure_HMA(ibox)%lattice = pressure(ibox)%computed - pressure(ibox)%ideal + END IF END DO WRITE(logunit,'(A80)') '********************************************************************************' diff --git a/Src/read_write_checkpoint.f90 b/Src/read_write_checkpoint.f90 index 16230d01..48864b60 100755 --- a/Src/read_write_checkpoint.f90 +++ b/Src/read_write_checkpoint.f90 @@ -441,6 +441,9 @@ SUBROUTINE Read_Config(ibox) atom_list(ia,this_im,is)%rxp, & atom_list(ia,this_im,is)%ryp, & atom_list(ia,this_im,is)%rzp + atom_list(ia,this_im,is)%rxp_init = atom_list(ia,this_im,is)%rxp + atom_list(ia,this_im,is)%ryp_init = atom_list(ia,this_im,is)%ryp + atom_list(ia,this_im,is)%rzp_init = atom_list(ia,this_im,is)%rzp ! set the frac and exist flags for this atom molecule_list(this_im,is)%frac = this_lambda atom_list(ia,this_im,is)%exist = .TRUE. diff --git a/Src/simulation_properties.f90 b/Src/simulation_properties.f90 index ca2bbe70..01f4356e 100755 --- a/Src/simulation_properties.f90 +++ b/Src/simulation_properties.f90 @@ -204,6 +204,59 @@ SUBROUTINE Compute_Pressure(this_box) END SUBROUTINE Compute_Pressure +! !**************************************************************************** + + SUBROUTINE Compute_HMA(this_box) + !*************************************************************************** + ! + ! This subroutine calculates the pressure of the the box + ! + ! CALLED BY: + ! Write_Properties + ! + ! CALLS : + ! Compute_System_Total_Force + ! + ! Written by Andrew Schultz on 02/10/21 + ! + !*************************************************************************** + + IMPLICIT NONE + + INTEGER :: this_box + + REAL(DP) :: dxp, dyp, dzp + REAL(DP) :: dx, dy, dz + Real(DP) :: fV + + INTEGER :: is, mol_is, im, ia + INTEGER :: this_nmols + + this_nmols = SUM(nmols(:,this_box)) + + energy_HMA(this_box)%harmonic = 1.5_DP / beta(this_box) & + * (this_nmols-1) + + ! compute Fdr + CALL Compute_System_Total_Force(this_box) + + energy_HMA(this_box)%anharmonic = 0.5_DP * energy_HMA(this_box)%sum_Fdr & + + energy(this_box)%total & + - energy_HMA(this_box)%lattice + energy_HMA(this_box)%total = 0.5_DP * energy_HMA(this_box)%sum_Fdr & + + energy(this_box)%total & + + energy_HMA(this_box)%harmonic + + fV = (beta(this_box)*pressure_HMA(this_box)%harmonic - this_nmols / box_list(this_box)%volume) & + / (3*(this_nmols)-1) + + pressure_HMA(this_box)%anharmonic = pressure(this_box)%computed - pressure(this_box)%ideal & + - pressure_HMA(this_box)%lattice + fV * energy_HMA(this_box)%sum_Fdr + pressure_HMA(this_box)%total = pressure_HMA(this_box)%harmonic + pressure(this_box)%computed & + - pressure(this_box)%ideal + fV * energy_HMA(this_box)%sum_Fdr + + END SUBROUTINE Compute_HMA + ! !**************************************************************************** END MODULE Simulation_Properties diff --git a/Src/type_definitions.f90 b/Src/type_definitions.f90 index fb140914..1f74401b 100755 --- a/Src/type_definitions.f90 +++ b/Src/type_definitions.f90 @@ -204,7 +204,7 @@ MODULE Type_Definitions ! atom_list has dimensions (natoms, max_molecules, nspecies) REAL(DP) :: rxp, ryp, rzp - REAL(DP) :: rxp_nls, ryp_nls, rzp_nls ! The starting positions for the neighbor list + REAL(DP) :: rxp_init, ryp_init, rzp_init ! The starting positions REAL(DP) :: rxp_old, ryp_old, rzp_old LOGICAL :: exist @@ -593,6 +593,53 @@ MODULE Type_Definitions END TYPE Pressure_Class +!------------------------------------------------------------------------------------------------- + + TYPE Energy_HMA_Class + ! This class holds the HMA energy values for each box + + ! Lattice energy -- energy at the beginning of the simulation + REAL(DP) :: lattice + + ! Harmonic energy: (3/2) (N-1) kT + REAL(DP) :: harmonic + + ! Anharmonic energy (from HMA) + REAL(DP) :: anharmonic + + ! Total energy (from HMA): U_lat + H_harm + U_anh + REAL(DP) :: total + + ! Sum of Fi dot dri over all atoms + REAL(DP) :: sum_Fdr + + ! last calculation + INTEGER :: last_calc + + END TYPE Energy_HMA_Class + +!------------------------------------------------------------------------------------------------- + + TYPE Pressure_HMA_Class + ! This class holds the HMA energy values for each box + + ! Lattice pressure -- non-ideal pressure at the beginning of the simulation + REAL(DP) :: lattice + + ! Harmonic pressure: read from input file + REAL(DP) :: harmonic + + ! Anharmonic pressure (from HMA) + REAL(DP) :: anharmonic + + ! Total pressure (from HMA): P_lat + P_harm + P_anh + REAL(DP) :: total + + ! last calculation + INTEGER :: last_calc + + END TYPE Pressure_HMA_Class + !------------------------------------------------------------------------------------------------- TYPE Rotation_Class diff --git a/Src/write_properties.f90 b/Src/write_properties.f90 index a8c51104..d79fad8f 100755 --- a/Src/write_properties.f90 +++ b/Src/write_properties.f90 @@ -148,6 +148,14 @@ SUBROUTINE Write_Header(file_number) prop_unit(ii) = '(kJ/mol)' + ELSE IF (prop_to_write(1:10) == 'Energy_HMA') THEN + + prop_unit(ii) = '(kJ/mol)-Ext' + + ELSE IF (prop_to_write(1:12) == 'Pressure_HMA') THEN + + prop_unit(ii) = '(bar)' + END IF WRITE(this_unit,'(A16,2X)',ADVANCE='NO') (TRIM(prop_unit(ii))) @@ -418,6 +426,29 @@ SUBROUTINE Write_Properties_Buffer(file_number) write_buff(ii+1) = mass_density / box_list(this_box)%volume END IF write_buff(ii+1) = write_buff(ii+1) * atomic_to_kgm3 + + ELSE IF (prop_written == 'Energy_HMA') THEN + IF (block_avg) THEN + write_buff(ii+1) = ac_energy_HMA(this_box,iblock) + ELSE + IF (energy_HMA(this_box)%last_calc /= i_mcstep) THEN + energy_HMA(this_box)%last_calc = i_mcstep + CALL Compute_HMA(this_box) + END IF + write_buff(ii+1) = energy_HMA(this_box)%total + END IF + write_buff(ii+1) = write_buff(ii+1) * atomic_to_kJmol + ELSE IF (prop_written == 'Pressure_HMA') THEN + IF (block_avg) THEN + write_buff(ii+1) = ac_pressure_HMA(this_box,iblock) + ELSE + IF (pressure_HMA(this_box)%last_calc /= i_mcstep) THEN + pressure_HMA(this_box)%last_calc = i_mcstep + CALL Compute_HMA(this_box) + END IF + write_buff(ii+1) = pressure_HMA(this_box)%total + END IF + write_buff(ii+1) = write_buff(ii+1) * atomic_to_bar END IF ! At the end increment property counter by 1 @@ -785,6 +816,24 @@ SUBROUTINE Write_Mean_Error(ibox) write_mean(ii+1) = write_mean(ii+1) * atomic_to_kgm3 write_err(ii+1) = sqrt(write_err(ii+1) / nbr_blocks) * atomic_to_kgm3 + ELSE IF (prop_written == 'Energy_HMA') THEN + + write_mean(ii+1) = SUM(ac_energy_HMA(ibox,:)) / nbr_blocks + DO iblock = 1, nbr_blocks + write_err(ii+1) = write_err(ii+1) + (ac_energy_HMA(ibox,iblock) - write_mean(ii+1))**2 + END DO + write_mean(ii+1) = write_mean(ii+1) * atomic_to_kJmol + write_err(ii+1) = sqrt(write_err(ii+1) / nbr_blocks) * atomic_to_kJmol + + ELSE IF (prop_written == 'Pressure_HMA') THEN + + write_mean(ii+1) = SUM(ac_pressure_HMA(ibox,:)) / nbr_blocks + DO iblock = 1, nbr_blocks + write_err(ii+1) = write_err(ii+1) + (ac_pressure_HMA(ibox,iblock) - write_mean(ii+1))**2 + END DO + write_mean(ii+1) = write_mean(ii+1) * atomic_to_bar + write_err(ii+1) = sqrt(write_err(ii+1) / nbr_blocks) * atomic_to_bar + END IF ! At the end increment property counter by 1