diff --git a/GUIDE.md b/GUIDE.md index c6ac74e..3c1af18 100644 --- a/GUIDE.md +++ b/GUIDE.md @@ -12,20 +12,44 @@ Most of the Fortran codes use a `namelist` to input a few parameters from standa This gives an easy way to specify default values in the program itself, and to use a keyword-based syntax to specify values different from the default ones at run-time. The input file, or string, should usually begin with `&nml` and end with `/`. -As a minimum, the program will expect to read `&nml /` (an empty list), but to -change the parameters, typical input might be +As a minimum, the program will expect to read `&nml /` (an empty list), +and on a Linux/Unix/bash system a typical command line would be ``` -&nml nblock=20, nstep=10000, dt=0.002 / +echo '&nml /' | ./mc_nvt_lj ``` -and the `key=value` pairs may be set out on different lines if you wish. +or similar. +To change the parameters, the namelist input might include the new values like this +``` +&nml nblock=20, nstep=10000, temperature=2.0 / +``` +Alternatively the namelist may be supplied in a file, +for instance `run.inp` containing +``` +&nml +nblock=20, +nstep=10000, +temperature=2.0 +/ +``` +and the program run like this +``` +./mc_nvt_lj < run.inp +``` +As indicated, the `key=value` pairs may be set out on different lines if you wish. ## Initial Configuration Simulation runs for bulk liquids require a starting configuration which can usually be prepared using the `initialize` program (built, by the default SConstruct file, in `build_initialize/`). The default parameters produce an FCC configuration of 256 atoms at reduced density ρ=0.75, writing out just the positions (for an MC program) to a file `cnf.inp`. -If the parameter `velocities=.true.` is supplied, then positions and velocities are -written to the file, corresponding to a reduced temperature _T_ = 1.0. +You would run the program as described above, +for instance on a Linux/Unix/bash system like this +``` +echo '&nml /' | ./initialize +``` +If the parameter `velocities=.true.` is supplied within the namelist, +then positions and velocities are written to the file, +corresponding to a reduced temperature _T_ = 1.0. These values of ρ and _T_ (see below) lie in the liquid region of the Lennard-Jones phase diagram. Non-default values may, of course, be supplied for this or other models. The `cnf.inp` file may then be copied to the directory in which the run is carried out. @@ -44,6 +68,32 @@ scales the velocities to change the kinetic energy per atom by a specified amoun and/or the positions (and the box length) to change the density by a specified amount. You may prefer to write your own program or script to perform these types of operation. +## Visualizing configurations +Our simulation configuration files have a very simple format: +the first line is the number of molecules, +the second line is the periodic box length (we always use cubic boxes), +or occasionally the bond length (for chain molecules which are simulated without a box), +and the third and subsequent lines each contain the coordinates of one atom or molecule. +The first three numbers on each of these lines are always the (x,y,z) position, +in simulation units (e.g. Lennard-Jones σ=1). +The rest of the numbers on each line contain velocities, orientation vectors etc., +as appropriate. + +This format is not compatible with most molecular visualization programs, +such as [JMOL](http://jmol.sourceforge.net/) +or [VMD](http://www.ks.uiuc.edu/Research/vmd/). +However, conversion into the basic [XYZ](https://en.wikipedia.org/wiki/XYZ_file_format) format +is easily achieved using a simple program, script, or editor. For example, +on most Linux/Unix systems, the `awk` language will be available: +``` +awk '(NR==1) {print} (NR==2) {print "Comment line"} (NR>2) {printf "%5s%15.6f%15.6f%15.6f\n", ("Ar"),($1*3.4),($2*3.4),($3*3.4)}' cnf.inp > cnf.xyz +``` +This produces a file which should be recognized as a set of Argon atoms with positions +(in Angstroms) appropriate to their van der Waals diameter. +This can be read into a molecular visualizer, +which will typically have various options for representing the atoms. +No attempt is made here to represent the periodicity of the system in the `cnf.xyz` file. + ## Lennard-Jones simulation programs A large number of the examples simulate the Lennard-Jones liquid. Before discussing these in detail, @@ -1197,8 +1247,8 @@ Larger values of _P_ give energies closer to the exact quantum mechanical canoni For this simple model, exact results can also be calculated for the finite values of _P_ used in the simulation -* KS Schweizer, RM Stratt, D Chandler, PG Wolynes, _J Chem Phys,_ __75,__ 1347 (1981). -* M Takahashi, M Imada, _J Phys Soc Japan,_ __53,__ 3765 (1984). +* KS Schweizer, RM Stratt, D Chandler, PG Wolynes, _J Chem Phys,_ __75,__ 1347 (1981), +* M Takahashi, M Imada, _J Phys Soc Japan,_ __53,__ 3765 (1984), and a routine to evaluate these is included in the example. No special techniques are used to accelerate the simulation; @@ -1220,3 +1270,72 @@ _P_ | _E_ (MC) | _E_ (exact) 6 | 0.4694(6) | 0.468708 7 | 0.4778(6) | 0.477941 8 | 0.4846(9) | 0.484244 + +The program `qmc_pi_lj` applies the path-integral method to the Lennard-Jones fluid. +The simplest, primitive, algorithm is used, +together with the crudest estimators for energy and pressure. +The program uses +single-particle Monte Carlo moves for the individual beads in the ring polymer, +along with translations of the centre-of-mass of each polymer. +As mentioned in the text, there are many improvements of all these aspects +of the algorithm, which are recommended for production work. + +The program takes in configuration files `cnf##.inp` where the `##` reflects +the polymer bead number, in the range 1 to _P_. +These files have the same format as the classical Lennard-Jones configurations. +They may be prepared in the same way as usual, +from the `initialize` program, from a run of `mc_nvt_lj` +at the appropriate density and temperature, +or from a run of `qmc_pi_lj` for a different value of _P_. +It does no harm if these starting configurations are simply duplicates of each other, +provided that a preliminary run is conducted to allow the polymers to equilibrate, +after which all the output files `cnf##.out` may be renamed to `cnf##.inp`. + +For testing, we compare with a set of simulations of neon, + +* M Neumann, M Zoppi, _Phys Rev E,_ __65,__ 031203 (2002), + +which are mainly based on an empirical pair potential, +but include selected results for Lennard-Jones for the case _P_=32. +The LJ parameters for neon are ε=36.8K, σ=0.2789nm, atomic mass _m_=20.18u, +and hence a reduced de Boer parameter λ=0.092σ, +which is the default value in the program. +We choose their lowest-temperature state point, +(_T_,ρ)=(25.8K,36.28nm-3)=(0.701087,0.787069) in reduced LJ units. +We use _N_=108 atoms, compared with Neumann and Zoppi's _N_=256, +and our runs are five times shorter (10 blocks of 10000 sweeps); +these differences should not be critical. +The maximum centre-of-mass displacement is 0.1σ. +Because the intra-polymer springs increase in strength with _P_, +the maximum displacement parameter for individual bead moves +is reduced from 0.06σ for _P_=2 +down to 0.02σ for _P_=32. +These choices give reasonable acceptance ratios for both kinds of move. +In the table below we report: +the rms radius _R_ of the ring polymers, +the average spring potential _E_(spring) per atom, +the kinetic energy _KE_ per atom, +total energy _E_ per atom, and pressure _p_ +(the last two including the standard long-range correction +for the applied cutoff _R_c=2.5σ). + +_P_ | _R_ | _E_(spring) | _KE_ | _E_ | _p_ +------ | ------ | ------ | ------ | ------ | ------ +1 † | 0.0 | 0.0 | 1.052 | -4.692(1) | -0.756(5) +2 | 0.04043(1) | 0.8963(7) | 1.2070(7) | -4.454(1) | -0.408(6) +4 | 0.04942(1) | 2.906(1) | 1.301(1) | -4.320(3) | -0.231(9) +8 | 0.05181(2) | 7.073(3) | 1.340(3) | -4.261(3) | -0.150(8) +16 | 0.05247(2) | 15.479(4) | 1.347(4) | -4.252(4) | -0.148(5) +32 | 0.05261(4) | 32.287(8) | 1.365(8) | -4.233(8) | -0.140(6) +32 ‡ | 0.053 | 32.3008 | 1.352 | -4.227 | -0.039 + +† For completeness the _P_=1 runs were performed using `mc_nvt_lj`. + +‡ These results are taken from Table I of Neumann and Zoppi (2002), +converted into LJ reduced units. + +A drawback of this state point is that the pressure is negative, +suggesting instability with respect to phase separation. +Nonetheless we have seen no evidence of crystalline structure +in the simulated configurations, +and assume that the liquid phase is metastable. diff --git a/python_examples/GUIDE.md b/python_examples/GUIDE.md index 28ce35f..dfe1bd1 100644 --- a/python_examples/GUIDE.md +++ b/python_examples/GUIDE.md @@ -94,11 +94,30 @@ but Python does not have this. Instead, to provide a keyword based syntax, we input these values using the widespread [JSON](http://www.json.org/ "JSON home page") format. -Typical input for a very simple example might be +As a minimum, the program will expect to read `{}` (an empty list), +and on a Linux/Unix/bash system a typical command line would be ``` -{ "nblock":10, "nstep":1000, "dt":0.005 } +echo '{}' | ./mc_nvt_lj.py ``` -and the `"key":value` pairs may be set out on different lines if you wish. +or similar. +To change the parameters, the JSON list might include the new values like this +``` +{ "nblock":20, "nstep":1000, "temperature":2.0 } +``` +Alternatively the list may be supplied in a file, +for instance `run.inp` containing +``` +{ +"nblock":20, +"nstep":10000, +"temperature":2.0 +} +``` +and the program run like this +``` +./mc_nvt_lj.py < run.inp +``` +As indicated, the `"key":value` pairs may be set out on different lines if you wish. The appearance is very similar to a Python dictionary, and indeed the data is loaded and parsed into a dictionary for further processing. Note carefully that we use a colon `:` rather than an equals sign to separate the @@ -122,8 +141,14 @@ Simulation runs for bulk liquids require a starting configuration which can usua the `initialize.py` program. The default parameters produce an FCC configuration of 256 atoms at reduced density ρ=0.75, writing out just the positions (for an MC program) to a file `cnf.inp`. -If the parameter `"velocities":true` is supplied, then positions and velocities are -written to the file, corresponding to a reduced temperature _T_ = 1.0. +You would run the program as described above, +for instance on a Linux/Unix/bash system like this +``` +echo '{}' | ./initialize.py +``` +If the parameter `"velocities":true` is supplied within the JSON list, +then positions and velocities are written to the file, +corresponding to a reduced temperature _T_ = 1.0. These values of ρ and _T_ (see below) lie in the liquid region of the Lennard-Jones phase diagram. Non-default values may, of course, be supplied for this or other models. The `cnf.inp` file may then be copied to the directory in which the run is carried out. @@ -143,6 +168,32 @@ scales the velocities to change the kinetic energy per atom by a specified amoun and/or the positions (and the box length) to change the density by a specified amount. You may prefer to write your own program or script to perform these types of operation. +## Visualizing configurations +Our simulation configuration files have a very simple format: +the first line is the number of molecules, +the second line is the periodic box length (we always use cubic boxes), +or occasionally the bond length (for chain molecules which are simulated without a box), +and the third and subsequent lines each contain the coordinates of one atom or molecule. +The first three numbers on each of these lines are always the (x,y,z) position, +in simulation units (e.g. Lennard-Jones σ=1). +The rest of the numbers on each line contain velocities, orientation vectors etc., +as appropriate. + +This format is not compatible with most molecular visualization programs, +such as [JMOL](http://jmol.sourceforge.net/) +or [VMD](http://www.ks.uiuc.edu/Research/vmd/). +However, conversion into the basic [XYZ](https://en.wikipedia.org/wiki/XYZ_file_format) format +is easily achieved using a simple program, script, or editor. For example, +on most Linux/Unix systems, the `awk` language will be available: +``` +awk '(NR==1) {print} (NR==2) {print "Comment line"} (NR>2) {printf "%5s%15.6f%15.6f%15.6f\n", ("Ar"),($1*3.4),($2*3.4),($3*3.4)}' cnf.inp > cnf.xyz +``` +This produces a file which should be recognized as a set of Argon atoms with positions +(in Angstroms) appropriate to their van der Waals diameter. +This can be read into a molecular visualizer, +which will typically have various options for representing the atoms. +No attempt is made here to represent the periodicity of the system in the `cnf.xyz` file. + ## Lennard-Jones simulation programs State points for the Lennard-Jones potential, in its cut (but not shifted) form, denoted (c), @@ -1059,3 +1110,75 @@ _P_ | _E_ (MC) | _E_ (exact) 6 | 0.4686(4) | 0.468708 7 | 0.4777(6) | 0.477941 8 | 0.4847(9) | 0.484244 + +The program `qmc_pi_lj.py` applies the path-integral method to the Lennard-Jones fluid. +The simplest, primitive, algorithm is used, +together with the crudest estimators for energy and pressure. +The program uses +single-particle Monte Carlo moves for the individual beads in the ring polymer, +along with translations of the centre-of-mass of each polymer. +As mentioned in the text, there are many improvements of all these aspects +of the algorithm, which are recommended for production work. + +The program takes in configuration files `cnf##.inp` where the `##` reflects +the polymer bead number, in the range 0 to _P_-1. +These files have the same format as the classical Lennard-Jones configurations. +They may be prepared in the same way as usual, +from the `initialize.py` program, from a run of `mc_nvt_lj.py` +at the appropriate density and temperature, +or from a run of `qmc_pi_lj.py` for a different value of _P_. +It does no harm if these starting configurations are simply duplicates of each other, +provided that a preliminary run is conducted to allow the polymers to equilibrate, +after which all the output files `cnf##.out` may be renamed to `cnf##.inp`. + +For testing, we compare with a set of simulations of neon, + +* M Neumann, M Zoppi, _Phys Rev E,_ __65,__ 031203 (2002), + +which are mainly based on an empirical pair potential, +but include selected results for Lennard-Jones for the case _P_=32. +The LJ parameters for neon are ε=36.8K, σ=0.2789nm, atomic mass _m_=20.18u, +and hence a reduced de Boer parameter λ=0.092σ, +which is the default value in the program. +We choose their lowest-temperature state point, +(_T_,ρ)=(25.8K,36.28nm-3)=(0.701087,0.787069) in reduced LJ units. +We use _N_=108 atoms, compared with Neumann and Zoppi's _N_=256, +and our runs are five times shorter (10 blocks of 10000 sweeps); +these differences should not be critical. +(Note that the program default run length is still shorter). +The maximum centre-of-mass displacement is 0.1σ. +Because the intra-polymer springs increase in strength with _P_, +the maximum displacement parameter for individual bead moves +is reduced from 0.06σ for _P_=2 +down to 0.03σ for _P_=16. +These choices give reasonable acceptance ratios for both kinds of move. +Because of the slowness of the program, we do not attempt the _P_=32 case; +nonetheless the trend is fairly clear in the table below, +and we can compare directly with the results of the corresponding Fortran example. +In the table below we report: +the rms radius _R_ of the ring polymers, +the average spring potential _E_(spring) per atom, +the kinetic energy _KE_ per atom, +total energy _E_ per atom, and pressure _p_ +(the last two including the standard long-range correction +for the applied cutoff _R_c=2.5σ). + +_P_ | _R_ | _E_(spring) | _KE_ | _E_ | _p_ +------ | ------ | ------ | ------ | ------ | ------ +1 † | 0.0 | 0.0 | 1.052 | -4.692(1) | -0.756(6) +2 | 0.04043(1) | 0.8961(6) | 1.2072(6) | -4.455(2) | -0.411(7) +4 | 0.04942(1) | 2.906(1) | 1.300(1) | -4.321(1) | -0.236(9) +8 | 0.05183(1) | 7.075(2) | 1.338(2) | -4.266(3) | -0.166(8) +16 | 0.05246(3) | 15.480(3) | 1.346(3) | -4.258(4) | -0.17(1) +32 ‡ | 0.053 | 32.3008 | 1.352 | -4.227 | -0.039 + +† For completeness the _P_=1 runs were performed using `mc_nvt_lj.py`. + +‡ These results are taken from Table I of Neumann and Zoppi (2002), +converted into LJ reduced units. + +A drawback of this state point is that the pressure is negative, +suggesting instability with respect to phase separation. +Nonetheless we have seen no evidence of crystalline structure +in the simulated configurations, +and assume that the liquid phase is metastable. diff --git a/python_examples/qmc_pi_lj.py b/python_examples/qmc_pi_lj.py new file mode 100755 index 0000000..0b4bec0 --- /dev/null +++ b/python_examples/qmc_pi_lj.py @@ -0,0 +1,293 @@ +#!/usr/bin/env python3 +# qmc_pi_lj.py + +#------------------------------------------------------------------------------------------------# +# This software was written in 2016/17 # +# by Michael P. Allen / # +# and Dominic J. Tildesley ("the authors"), # +# to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"), # +# published by Oxford University Press ("the publishers"). # +# # +# LICENCE # +# Creative Commons CC0 Public Domain Dedication. # +# To the extent possible under law, the authors have dedicated all copyright and related # +# and neighboring rights to this software to the PUBLIC domain worldwide. # +# This software is distributed without any warranty. # +# You should have received a copy of the CC0 Public Domain Dedication along with this software. # +# If not, see . # +# # +# DISCLAIMER # +# The authors and publishers make no warranties about the software, and disclaim liability # +# for all uses of the software, to the fullest extent permitted by applicable law. # +# The authors and publishers do not recommend use of this software for any purpose. # +# It is made freely available, solely to clarify points made in the text. When using or citing # +# the software, you should not imply endorsement by the authors or publishers. # +#------------------------------------------------------------------------------------------------# + +"""Quantum Monte Carlo, path-integral method.""" + +def calc_variables ( ): + """Calculates all variables of interest. + + They are collected and returned as a list, for use in the main program. + """ + + # In this example we simulate using the cut (but not shifted) potential + # but we only report results which have had the long-range corrections applied + # The value of the cut-and-shifted potential is not used, in this example + + import numpy as np + import math + from averages_module import VariableType + from lrc_module import potential_lrc, pressure_lrc + + # Preliminary calculations (n,r,total are taken from the calling program) + vol = box**3 # Volume + rho = n / vol # Density + kin = 1.5 * n * p * temperature # Average kinetic energy for NP-atom system + kin_q = kin - total_spr # Quantum estimator for kinetic energy + rad_g = rad_gyr ( r ) + + # Variables of interest, of class VariableType, containing three attributes: + # .val: the instantaneous value + # .nam: used for headings + # .method: indicating averaging method + # If not set below, .method adopts its default value of avg + # The .nam and some other attributes need only be defined once, at the start of the program, + # but for clarity and readability we assign all the values together below + + # Acceptance ratio of atomic moves + r_r = VariableType ( nam = 'Atomic move ratio', val = r_ratio, instant = False ) + + # Acceptance ratio of centre-of-mass moves + c_r = VariableType ( nam = 'COM move ratio', val = c_ratio, instant = False ) + + # Internal energy per atom for full potential with LRC + # LRC plus cut (but not shifted) PE already divided by factor P + # plus KE estimator: total classical KE for NP-atom system MINUS total spring potential + # all divided by N + e_f = VariableType ( nam = 'E/N full', val = potential_lrc(rho,r_cut) + (kin_q+total.pot)/n ) + + # Kinetic energy per atom, just for interest's sake + k_q = VariableType ( nam = 'KE/N', val = kin_q/n ) + + # Pressure for full potential with LRC + # LRC plus ideal gas contribution plus total virial divided by V + kin_q = kin_q / 1.5 # Convert KE estimator to kinetic energy part of pressure + p_f = VariableType ( nam = 'P full', val = pressure_lrc(rho,r_cut) + (kin_q+total.vir)/vol ) + + # Quantum spring energy per atom, just for interest's sake + e_q = VariableType ( nam = 'Espring/N', val = total_spr/n ) + + # Quantum polymer radius of gyration, just for interest's sake + r_g = VariableType ( nam = 'Radius of gyration', val = rad_g ) + + # Collect together into a list for averaging + return [ r_r, c_r, e_f, p_f, e_q, k_q, r_g ] + +def rad_gyr ( r ): + """Calculate average radius of gyration of polymers.""" + + import numpy as np + + # The formula we use involves a single sweep over atoms, and is origin-independent + # To implement periodic boundaries, we take the origin on atom 1 + + r_g = 0.0 # Zero function accumulator + p, n, d = r.shape + assert d==3, 'Dimension error for r' + + r_g = 0.0 + + for i in range(n): # Loop over polymers + ri = r[1:,i,:] - r[0,i,:] # Get coordinates of beads 1.. in polymer relative to bead 0 + ri = ri - np.rint(ri) # Position with PBC applied (box = 1 units) + r_cm = np.sum(ri,axis=0) / p # centre-of-mass vector + r_sq = np.sum(ri**2) / p # mean-squared distance + r_sq = r_sq - np.sum(r_cm**2) # squared radius of gyration + if r_sq<0.0: + r_sq=0.0 # Guard against roundoff + r_g = r_g + np.sqrt(r_sq) # Accumulate root-mean-square radius of gyration + + return box * r_g / n # Average RMS Rg in sigma=1 units + +# Takes in a configuration of atoms (positions) +# Cubic periodic boundary conditions +# Conducts path-integral Monte Carlo at the given temperature +# Uses no special neighbour lists + +# Reads several variables and options from standard input using JSON format +# Leave input empty "{}" to accept supplied defaults + +# Positions r are divided by box length after reading in +# However, input configuration, output configuration, most calculations, and all results +# are given in simulation units defined by the model +# For example, for Lennard-Jones, sigma = 1, epsilon = 1 +# The importance of quantum effects is specified through the reduced de Boer length +# deboer = (hbar/sqrt(mass*epsilon))/sigma which takes typical values of +# 0.01 for Xe, 0.03 for Ar, and 0.095 for Ne. +# This means that the quantum spring constant may be expressed k_spring = P*(T/deboer)**2 +# where T stands for the reduced temperature kB*T/epsilon + +# Despite the program name, there is nothing here specific to Lennard-Jones +# The model is defined in qmc_module + +import json +import sys +import numpy as np +from config_io_module import read_cnf_atoms, write_cnf_atoms +from averages_module import run_begin, run_end, blk_begin, blk_end, blk_add +from maths_module import random_translate_vector, metropolis +from qmc_pi_lj_module import introduction, conclusion, potential, spring, potential_1, spring_1, PotentialType + +cnf_prefix = 'cnf' +inp_tag = 'inp' +out_tag = 'out' +sav_tag = 'sav' + +print('qmc_pi_lj') +print('Path-integral Monte Carlo, constant-NVT ensemble') +print('Simulation uses cut (but not shifted) potential') + +# Read parameters in JSON format +try: + nml = json.load(sys.stdin) +except json.JSONDecodeError: + print('Exiting on Invalid JSON format') + sys.exit() + +# Set default values, check keys and typecheck values +defaults = {"nblock":10, "nstep":1000, "temperature":0.701087, "r_cut":2.5, + "dr_max":0.05, "dc_max":0.1, "p":4, "deboer":0.092} +for key, val in nml.items(): + if key in defaults: + assert type(val) == type(defaults[key]), key+" has the wrong type" + else: + print('Warning', key, 'not in ',list(defaults.keys())) + +# Set parameters to input values or defaults +nblock = nml["nblock"] if "nblock" in nml else defaults["nblock"] +nstep = nml["nstep"] if "nstep" in nml else defaults["nstep"] +temperature = nml["temperature"] if "temperature" in nml else defaults["temperature"] +r_cut = nml["r_cut"] if "r_cut" in nml else defaults["r_cut"] +dr_max = nml["dr_max"] if "dr_max" in nml else defaults["dr_max"] +dc_max = nml["dc_max"] if "dc_max" in nml else defaults["dc_max"] +p = nml["p"] if "p" in nml else defaults["p"] +deboer = nml["deboer"] if "deboer" in nml else defaults["deboer"] + +introduction() +np.random.seed() + +# Write out parameters +print( "{:40}{:15d} ".format('Size of each ring polymer', p) ) +print( "{:40}{:15d} ".format('Number of blocks', nblock) ) +print( "{:40}{:15d} ".format('Number of steps per block', nstep) ) +print( "{:40}{:15.6f}".format('Temperature', temperature) ) +print( "{:40}{:15.6f}".format('Potential cutoff distance', r_cut) ) +print( "{:40}{:15.6f}".format('Maximum displacement', dr_max) ) +print( "{:40}{:15.6f}".format('Maximum COM displacement', dc_max) ) +print( "{:40}{:15.6f}".format('de Boer length', deboer) ) +assert p>1 and p<100, 'p must lie between 2 and 99' +k_spring = p * ( temperature / deboer ) ** 2 +print( "{:40}{:15.6f}".format('Quantum spring constant', k_spring) ) + +# Read in initial configuration +# Read each polymer index from a unique file; we assume that n and box are all the same! +n, box, r = read_cnf_atoms ( cnf_prefix+'00.'+inp_tag) # Read first LJ configuration to get basic parameters n and box +print( "{:40}{:15d} ".format('Number of particles', n) ) +print( "{:40}{:15.6f}".format('Box length', box) ) +print( "{:40}{:15.6f}".format('Density', n/box**3) ) +r = np.empty( (p,n,3), dtype=np.float_ ) +for k in range(p): + n, box, r[k,:,:] = read_cnf_atoms ( cnf_prefix+str(k).zfill(2)+'.'+inp_tag) +r = r / box # Convert positions to box units +r = r - np.rint ( r ) # Periodic boundaries + +# Calculate classical LJ and quantum spring potential energies & check overlap +total = potential ( box, r_cut, r ) +assert not total.ovr, 'Overlap in initial configuration' +total_spr = spring ( box, k_spring, r ) + +# Initialize arrays for averaging and write column headings +r_ratio = 0.0 +c_ratio = 0.0 +run_begin ( calc_variables() ) + +for blk in range(1,nblock+1): # Loop over blocks + + blk_begin() + + for stp in range(nstep): # Loop over steps + + r_moves = 0 + c_moves = 0 + + for i in range(n): # Loop over atoms + + # Centre of mass move + dc = random_translate_vector ( dc_max/box, np.zeros(3,dtype=np.float_) ) + partial_old = PotentialType ( 0.0, 0.0, False ) + partial_new = PotentialType ( 0.0, 0.0, False ) + for k in range(p): # Loop over ring polymer indices + rkj = np.delete(r[k,:,:],i,0) # Array of all other atoms for this k + partial_old = partial_old + potential_1 ( r[k,i,:], box, r_cut, rkj, p ) # Old atom classical potential etc + rki = r[k,i,:] + dc + rki = rki - np.rint(rki) + partial_new = partial_new + potential_1 ( rki, box, r_cut, rkj, p ) # New atom classical potential etc + assert not partial_old.ovr, 'Overlap in current configuration' + if not partial_new.ovr: # Test for non-overlapping configuration + delta = partial_new.pot - partial_old.pot # Change in classical cut (but not shifted) potential + delta = delta / temperature + if metropolis ( delta ): # Accept Metropolis test + total = total + partial_new - partial_old # Update total values + r[:,i,:] = r[:,i,:] + dc # Update positions + r[:,i,:] = r[:,i,:] - np.rint(r[:,i,:]) # Periodic boundary conditions + c_moves = c_moves + 1 # Increment move counter + + # Individual atom moves + for k in range(p): # Loop over ring polymer indices + rkj = np.delete(r[k,:,:],i,0) # Array of all the other atoms + partial_old = potential_1 ( r[k,i,:], box, r_cut, rkj, p ) # Old atom classical potential etc + assert not partial_old.ovr, 'Overlap in current configuration' + + kp = (k+1)%p + km = (k-1)%p + partial_old_spr = ( spring_1 ( r[k,i,:], r[km,i,:], box, k_spring ) + + spring_1 ( r[k,i,:], r[kp,i,:], box, k_spring ) ) + + rki = random_translate_vector ( dr_max/box, r[k,i,:] ) # Trial move to new position (in box=1 units) + rki = rki - np.rint ( rki ) # Periodic boundary correction + partial_new = potential_1 ( rki, box, r_cut, rkj, p ) # New atom potential, virial etc + + if not partial_new.ovr: # Test for non-overlapping configuration + partial_new_spr = ( spring_1 ( rki, r[km,i,:], box, k_spring ) + + spring_1 ( rki, r[kp,i,:], box, k_spring ) ) + delta = partial_new.pot - partial_old.pot # Change in classical cut (but not shifted) potential + delta = delta + partial_new_spr - partial_old_spr # Add change in quantum potential + delta = delta / temperature + if metropolis ( delta ): # Accept Metropolis test + total = total + partial_new - partial_old # Update total values + total_spr = total_spr + partial_new_spr - partial_old_spr # Update quantum system potential + r[k,i,:] = rki # Update position + r_moves = r_moves + 1 # Increment move counter + + r_ratio = r_moves / (n*p) + c_ratio = c_moves / n + + blk_add ( calc_variables() ) + + blk_end(blk) # Output block averages + sav_tag = str(blk).zfill(3) if blk<1000 else 'sav' # Number configuration by block + for k in range(p): + write_cnf_atoms ( cnf_prefix+str(k).zfill(2)+'.'+sav_tag, n, box, r[k,:,:]*box ) # Save configuration + +run_end ( calc_variables() ) + +total = potential ( box, r_cut, r ) # Double check book-keeping +assert not total.ovr, 'Overlap in final configuration' +total_spr = spring ( box, k_spring, r ) + +for k in range(p): + write_cnf_atoms ( cnf_prefix+str(k).zfill(2)+'.'+out_tag, n, box, r[k,:,:]*box ) # Save configuration + +conclusion() diff --git a/python_examples/qmc_pi_lj_module.py b/python_examples/qmc_pi_lj_module.py new file mode 100644 index 0000000..6b05a50 --- /dev/null +++ b/python_examples/qmc_pi_lj_module.py @@ -0,0 +1,192 @@ +#!/usr/bin/env python3 +# qmc_pi_lj_module.py + +#------------------------------------------------------------------------------------------------# +# This software was written in 2016/17 # +# by Michael P. Allen / # +# and Dominic J. Tildesley ("the authors"), # +# to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"), # +# published by Oxford University Press ("the publishers"). # +# # +# LICENCE # +# Creative Commons CC0 Public Domain Dedication. # +# To the extent possible under law, the authors have dedicated all copyright and related # +# and neighboring rights to this software to the PUBLIC domain worldwide. # +# This software is distributed without any warranty. # +# You should have received a copy of the CC0 Public Domain Dedication along with this software. # +# If not, see . # +# # +# DISCLAIMER # +# The authors and publishers make no warranties about the software, and disclaim liability # +# for all uses of the software, to the fullest extent permitted by applicable law. # +# The authors and publishers do not recommend use of this software for any purpose. # +# It is made freely available, solely to clarify points made in the text. When using or citing # +# the software, you should not imply endorsement by the authors or publishers. # +#------------------------------------------------------------------------------------------------# + +"""Energy and move routines for PIMC simulation, LJ potential.""" + +fast = True # Change this to replace NumPy potential evaluation with slower Python + +class PotentialType: + """A composite variable for interactions.""" + + def __init__(self, pot, vir, ovr): + self.pot = pot # the potential energy cut (but not shifted) at r_cut + self.vir = vir # the virial + self.ovr = ovr # a flag indicating overlap (i.e. pot too high to use) + + def __add__(self, other): + pot = self.pot + other.pot + vir = self.vir + other.vir + ovr = self.ovr or other.ovr + return PotentialType(pot,vir,ovr) + + def __sub__(self, other): + pot = self.pot - other.pot + vir = self.vir - other.vir + ovr = self.ovr or other.ovr # This is meaningless, but inconsequential + return PotentialType(pot,vir,ovr) + +def introduction(): + """Prints out introductory statements at start of run.""" + + print('Lennard-Jones potential') + print('Cut (but not shifted)') + print('Diameter, sigma = 1') + print('Well depth, epsilon = 1') + if fast: + print('Fast NumPy potential routine') + else: + print('Slow Python potential routine') + +def conclusion(): + """Prints out concluding statements at end of run.""" + + print('Program ends') + +def potential ( box, r_cut, r ): + """Takes in box, cutoff range, and coordinate array, and calculates total potential etc. + + The results are returned as total, a PotentialType variable. + """ + # Actual calculation performed by function potential_1 + + p, n, d = r.shape + assert d==3, 'Dimension error for r' + + total = PotentialType ( pot=0.0, vir=0.0, ovr=False ) + + for k in range(p): + for i in range(n-1): + partial = potential_1 ( r[k,i,:], box, r_cut, r[k,i+1:,:], p ) + if partial.ovr: + total.ovr = True + return total + total = total + partial + + return total + +def potential_1 ( rki, box, r_cut, r, p ): + """Takes in coordinates of an atom and calculates its interactions. + + Values of box, cutoff range, and partner coordinate array are supplied. + The results are returned as partial, a PotentialType variable. + """ + + import numpy as np + + # partial.pot is the nonbonded cut (not shifted) potential energy of atom rki with a set of other atoms + # partial.vir is the corresponding virial of atom rki + # partial.lap is the corresponding Laplacian of atom rki + # partial.ovr is a flag indicating overlap (potential too high) to avoid overflow + # If this is True, the values of partial.pot etc should not be used + # In general, r will be a subset of the complete set of simulation coordinates + # and none of its rows should be identical to rki + + # It is assumed that positions are in units where box = 1 + # Forces are calculated in units where sigma = 1 and epsilon = 1 + + nj, d = r.shape + assert d==3, 'Dimension error for r in potential_1' + assert rki.size==3, 'Dimension error for rki in potential_1' + + sr2_ovr = 1.77 # Overlap threshold (pot > 100) + r_cut_box = r_cut / box + r_cut_box_sq = r_cut_box ** 2 + box_sq = box ** 2 + + if fast: + r_ki_kj = rki - r # Get all separation vectors from partners + r_ki_kj = r_ki_kj - np.rint(r_ki_kj) # Periodic boundary conditions in box=1 units + r_ki_kj_sq = np.sum(r_ki_kj**2,axis=1) # Squared separations + in_range = r_ki_kj_sq < r_cut_box_sq # Set flags for within cutoff + r_ki_kj_sq = r_ki_kj_sq * box_sq # Now in sigma=1 units + sr2 = np.where ( in_range, 1.0 / r_ki_kj_sq, 0.0 ) # (sigma/r_ki_kj)**2, only if in range + ovr = sr2 > sr2_ovr # Set flags for any overlaps + if np.any(ovr): + partial = PotentialType ( pot=0.0, vir=0.0, ovr=True ) + return partial + sr6 = sr2 ** 3 + sr12 = sr6 ** 2 + pot = sr12 - sr6 # LJ pair potentials (cut but not shifted) + vir = pot + sr12 # LJ pair virials + partial = PotentialType ( pot=np.sum(pot), vir=np.sum(vir), ovr=False ) + + else: + partial = PotentialType ( pot=0.0, vir=0.0, ovr=False ) + for rkj in r: + r_ki_kj = rki - rkj # Separation vector + r_ki_kj = r_ki_kj - np.rint(r_ki_kj) # Periodic boundary conditions in box=1 units + r_ki_kj_sq = np.sum(r_ki_kj**2) # Squared separation + if r_ki_kj_sq < r_cut_box_sq: # Check within cutoff + r_ki_kj_sq = r_ki_kj_sq * box_sq # Now in sigma=1 units + sr2 = 1.0 / r_ki_kj_sq # (sigma/r_ki_kj)**2 + ovr = sr2 > sr2_ovr # Overlap if too close + if ovr: + partial.ovr=True + return partial + sr6 = sr2 ** 3 + sr12 = sr6 ** 2 + pot = sr12 - sr6 # LJ pair potential (cut but not shifted) + vir = pot + sr12 # LJ pair virial + partial = partial + PotentialType ( pot=pot, vir=vir, ovr=ovr ) + + # Include numerical factors + partial.pot = partial.pot * 4.0 / p # Classical potentials are weaker by a factor p + partial.vir = partial.vir * 24.0 / (3*p) # 24*epsilon and divide virial by 3 & by p + + return partial + +def spring ( box, k_spring, r ): + """Takes in box, spring strength, and coordinate array, and calculates quantum spring potential.""" + + import numpy as np + + # Actual calculation performed by function spring_1 + + p, n, d = r.shape + assert d==3, 'Dimension error for r in potential' + + total = 0.0 + + for i in range(n): # Loop over ring polymers + for k in range(p): # Loop over atoms within polymer + kp = (k+1)%p # Neighbour index + partial = spring_1 ( r[k,i,:], r[kp,i,:], box, k_spring ) + total = total + partial + + return total + +def spring_1 ( rki, rli, box, k_spring ): + """Returns quantum potential for given atom.""" + + import numpy as np + + # Coordinate array of neighbour is supplied + + r_ki_li = rki - rli # Separation vector + r_ki_li = r_ki_li - np.rint(r_ki_li) # Periodic boundary conditions in box=1 units + r_ki_li_sq = np.sum(r_ki_li**2) * box**2 # Squared separation in sigma=1 units + + return 0.5 * k_spring * r_ki_li_sq diff --git a/qmc_pi_lj.f90 b/qmc_pi_lj.f90 index e56bb2b..a56485e 100644 --- a/qmc_pi_lj.f90 +++ b/qmc_pi_lj.f90 @@ -59,15 +59,18 @@ PROGRAM qmc_pi_lj REAL :: box ! Box length REAL :: lambda ! de Boer length REAL :: dr_max ! Maximum MC displacement + REAL :: dc_max ! Maximum centre-of-mass MC displacement REAL :: temperature ! Specified temperature REAL :: r_cut ! Potential cutoff distance ! Composite interaction = pot & ovr variables TYPE(potential_type) :: total, partial_old, partial_new - INTEGER :: blk, stp, i, k, nstep, nblock, moves, ioerr - REAL :: total_spr, partial_old_spr, partial_new_spr, delta, k_spring, m_ratio - REAL, DIMENSION(3) :: rik + INTEGER :: blk, stp, i, k, km, kp, nstep, nblock, r_moves, c_moves, ioerr + REAL :: total_spr, partial_old_spr, partial_new_spr, delta, k_spring + REAL :: r_ratio, c_ratio + REAL :: rad_g + REAL, DIMENSION(3) :: rik, dc CHARACTER(len=3), PARAMETER :: inp_tag = 'inp' CHARACTER(len=3), PARAMETER :: out_tag = 'out' @@ -75,7 +78,7 @@ PROGRAM qmc_pi_lj CHARACTER(len=6) :: cnf_prefix = 'cnf##.' ! Will have polymer id inserted CHARACTER(len=2) :: k_tag ! Will contain polymer id - NAMELIST /nml/ nblock, nstep, p, temperature, r_cut, dr_max, lambda + NAMELIST /nml/ nblock, nstep, p, temperature, r_cut, dr_max, dc_max, lambda WRITE ( unit=output_unit, fmt='(a)' ) 'qmc_pi_lj' WRITE ( unit=output_unit, fmt='(a)' ) 'Path-integral Monte Carlo, constant-NVT ensemble' @@ -86,12 +89,13 @@ PROGRAM qmc_pi_lj ! Set sensible default run parameters for testing nblock = 10 - nstep = 1000 - temperature = 0.7 + nstep = 10000 + temperature = 0.701087 r_cut = 2.5 - dr_max = 0.15 + dr_max = 0.05 + dc_max = 0.1 p = 4 - lambda = 0.1 + lambda = 0.092 ! from Ne parameters epsilon=36.8K, sigma=0.2789nm, m=20.18u ! Read run parameters from namelist ! Comment out, or replace, this section if you don't like namelists @@ -110,6 +114,7 @@ PROGRAM qmc_pi_lj WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Temperature', temperature WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Potential cutoff distance', r_cut WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Maximum displacement', dr_max + WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Maximum COM displacement', dc_max WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'de Boer length', lambda IF ( p < 2 .OR. p > 99 ) THEN WRITE ( unit=error_unit, fmt='(a,i15)') 'p must lie between 2 and 99', p @@ -143,7 +148,8 @@ PROGRAM qmc_pi_lj total_spr = spring ( box, k_spring ) ! Initialize arrays for averaging and write column headings - m_ratio = 0.0 + r_ratio = 0.0 + c_ratio = 0.0 CALL run_begin ( calc_variables() ) DO blk = 1, nblock ! Begin loop over blocks @@ -152,10 +158,45 @@ PROGRAM qmc_pi_lj DO stp = 1, nstep ! Begin loop over steps - moves = 0 + r_moves = 0 + c_moves = 0 DO i = 1, n ! Begin loop over atoms + ! Centre of mass move + + dc = 0.0 + dc = random_translate_vector ( dc_max/box, dc ) ! Random displacement + partial_old = potential_type ( 0.0, 0.0, .FALSE. ) + partial_new = potential_type ( 0.0, 0.0, .FALSE. ) + DO k = 1, p ! Loop over ring polymer indices + partial_old = partial_old + potential_1 ( r(:,i,k), i, k, box, r_cut ) ! Old atom classical potential etc + rik(:) = r(:,i,k) + dc + rik(:) = rik(:) - ANINT ( rik(:) ) ! Periodic boundary correction + partial_new = partial_new + potential_1 ( rik, i, k, box, r_cut ) ! New atom classical potential etc + END DO + + IF ( partial_old%ovr ) THEN ! should never happen + WRITE ( unit=error_unit, fmt='(a)') 'Overlap in current configuration' + STOP 'Error in qmc_pi_lj' + END IF + + IF ( .NOT. partial_new%ovr ) THEN ! Test for non-overlapping configuration + + delta = partial_new%pot - partial_old%pot ! Change in classical cut (but not shifted) potential + delta = delta / temperature ! Divide by temperature + + IF ( metropolis ( delta ) ) THEN ! Accept Metropolis test + total = total + partial_new - partial_old ! Update classical total values + r(:,i,:) = r(:,i,:) + SPREAD ( dc, dim=2, ncopies=p ) ! Update position + r(:,i,:) = r(:,i,:) - ANINT ( r(:,i,:) ) ! Periodic boundary corrections + c_moves = c_moves + 1 ! Increment move counter + END IF ! End accept Metropolis test + + END IF ! End test for non-overlapping configuration + + ! Individual atom moves + DO k = 1, p ! Begin loop over ring polymer indices partial_old = potential_1 ( r(:,i,k), i, k, box, r_cut ) ! Old atom classical potential etc @@ -165,7 +206,13 @@ PROGRAM qmc_pi_lj STOP 'Error in qmc_pi_lj' END IF - partial_old_spr = spring_1 ( r(:,i,k), i, k, box, k_spring ) ! Old atom quantum potential + ! Get neighbour indices + kp = k+1 + IF ( kp > p ) kp = 1 + km = k-1 + if ( km < 1 ) km = p + partial_old_spr = spring_1 ( r(:,i,k), r(:,i,km), box, k_spring ) & + & + spring_1 ( r(:,i,k), r(:,i,kp), box, k_spring ) ! Old atom quantum potential rik(:) = random_translate_vector ( dr_max/box, r(:,i,k) ) ! Trial move to new position (in box=1 units) rik(:) = rik(:) - ANINT ( rik(:) ) ! Periodic boundary correction @@ -174,7 +221,8 @@ PROGRAM qmc_pi_lj IF ( .NOT. partial_new%ovr ) THEN ! Test for non-overlapping configuration - partial_new_spr = spring_1 ( rik, i, k, box, k_spring ) ! New atom quantum potential + partial_new_spr = spring_1 ( rik, r(:,i,km), box, k_spring ) & + & + spring_1 ( rik, r(:,i,kp), box, k_spring ) ! New atom quantum potential delta = partial_new%pot - partial_old%pot ! Change in classical cut (but not shifted) potential delta = delta + partial_new_spr - partial_old_spr ! Add change in quantum potential @@ -184,7 +232,7 @@ PROGRAM qmc_pi_lj total = total + partial_new - partial_old ! Update classical total values total_spr = total_spr + partial_new_spr - partial_old_spr ! Update quantum system potential r(:,i,k) = rik ! Update position - moves = moves + 1 ! Increment move counter + r_moves = r_moves + 1 ! Increment move counter END IF ! End accept Metropolis test END IF ! End test for non-overlapping configuration @@ -193,7 +241,8 @@ PROGRAM qmc_pi_lj END DO ! End loop over atoms - m_ratio = REAL(moves) / REAL(n) + r_ratio = REAL(r_moves) / REAL(n*p) + c_ratio = REAL(c_moves) / REAL(n) ! Calculate and accumulate variables for this step CALL blk_add ( calc_variables() ) @@ -232,27 +281,27 @@ PROGRAM qmc_pi_lj CONTAINS FUNCTION calc_variables ( ) RESULT ( variables ) - USE lrc_module, ONLY : potential_lrc + USE lrc_module, ONLY : potential_lrc, pressure_lrc USE averages_module, ONLY : variable_type IMPLICIT NONE - TYPE(variable_type), DIMENSION(3) :: variables ! The 3 variables listed below + TYPE(variable_type), DIMENSION(7) :: variables ! The 7 variables listed below ! This routine calculates all variables of interest and (optionally) writes them out ! They are collected together in the variables array, for use in the main program ! In this example we simulate using the cut (but not shifted) potential - ! The values of < e_c > and < density > should be consistent (for this potential) - ! For comparison, long-range corrections are also applied to give - ! an estimate of < e_f > for the full (uncut) potential + ! but we only report results which have had the long-range corrections applied ! The value of the cut-and-shifted potential is not used, in this example - TYPE(variable_type) :: m_r, e_c, e_f - REAL :: vol, rho, kin + TYPE(variable_type) :: r_r, c_r, e_f, p_f, e_q, k_q, r_g + REAL :: vol, rho, kin, kin_q ! Preliminary calculations - vol = box**3 ! Volume - rho = REAL(n) / vol ! Density - kin = 1.5 * n * p * temperature ! Average kinetic energy for NP-atom system + vol = box**3 ! Volume + rho = REAL(n) / vol ! Density + kin = 1.5 * n * p * temperature ! Average kinetic energy for NP-atom system + kin_q = kin - total_spr ! Quantum estimator for kinetic energy + rad_g = rad_gyr ( r ) ! Variables of interest, of type variable_type, containing three components: ! %val: the instantaneous value @@ -262,25 +311,74 @@ FUNCTION calc_variables ( ) RESULT ( variables ) ! The %nam and some other components need only be defined once, at the start of the program, ! but for clarity and readability we assign all the values together below - ! Acceptance ratio of moves - m_r = variable_type ( nam = 'Move ratio', val = m_ratio, instant = .FALSE. ) + ! Acceptance ratio of atomic moves + r_r = variable_type ( nam = 'Atomic move ratio', val = r_ratio, instant = .FALSE. ) - ! Internal energy per atom for simulated, cut, potential - ! Total (cut but not shifted) PE already divided by factor P - ! plus total classical KE for NP-atom system MINUS total spring potential - ! all divided by N - e_c = variable_type ( nam = 'E/N cut', val = (kin+total%pot-total_spr)/REAL(n) ) + ! Acceptance ratio of centre-of-mass moves + c_r = variable_type ( nam = 'COM move ratio', val = c_ratio, instant = .FALSE. ) ! Internal energy per atom for full potential with LRC ! LRC plus total (cut but not shifted) PE already divided by factor P - ! plus total classical KE for NP-atom system MINUS total spring potential + ! plus KE estimator: total classical KE for NP-atom system MINUS total spring potential ! all divided by N - e_f = variable_type ( nam = 'E/N full', val = potential_lrc(rho,r_cut) + (kin+total%pot-total_spr)/REAL(n) ) + e_f = variable_type ( nam = 'E/N full', val = potential_lrc(rho,r_cut) + (kin_q+total%pot)/REAL(n) ) + + ! Kinetic energy per atom, just for interest's sake + k_q = variable_type ( nam = 'KE/N', val = kin_q/REAL(n) ) + + ! Pressure for full potential with LRC + ! LRC plus ideal gas contribution plus total virial divided by V + kin_q = kin_q / 1.5 ! Convert KE estimator to kinetic energy part of pressure + p_f = variable_type ( nam = 'P full', val = pressure_lrc(rho,r_cut) + (kin_q+total%vir)/vol ) + + ! Quantum spring energy per atom, just for interest's sake + e_q = variable_type ( nam = 'Espring/N', val = total_spr/REAL(n) ) + + ! Quantum polymer radius of gyration, just for interest's sake + r_g = variable_type ( nam = 'Radius of gyration', val = rad_g ) ! Collect together for averaging - variables = [ m_r, e_c, e_f ] + variables = [ r_r, c_r, e_f, p_f, e_q, k_q, r_g ] END FUNCTION calc_variables + FUNCTION rad_gyr ( r ) RESULT ( r_g ) + IMPLICIT NONE + REAL, DIMENSION(:,:,:), INTENT(in) :: r ! Coordinate array + REAL :: r_g ! RMS radius of gyration + + ! Calculate average radius of gyration of polymers + ! The formula we use involves a single sweep over atoms, and is origin-independent + ! To implement periodic boundaries, we take the origin on atom 1 + + INTEGER :: i, k + REAL, DIMENSION(3) :: r_ik, r_cm + REAL :: r_sq + + r_g = 0.0 ! Zero function accumulator + + DO i = 1, n ! Loop over polymers + r_cm = 0.0 + r_sq = 0.0 + DO k = 2, p ! Loop over atoms in polymer, omitting first one + r_ik = r(:,i,k) - r(:,i,1) ! Position relative to atom 1 + r_ik = r_ik - ANINT(r_ik) ! Position with PBC applied (box = 1 units) + r_cm = r_cm + r_ik ! Increment centre-of-mass accumulator + r_sq = r_sq + SUM(r_ik**2) ! Increment squared distance accumulator + END DO ! End loop over atoms in polymer, omitting first one + r_cm = r_cm / REAL(p) ! Normalize centre-of-mass vector + r_sq = r_sq / REAL(p) ! Normalize mean-squared distances + r_sq = r_sq - SUM(r_cm**2) ! Definition of squared radius of gyration + IF ( r_sq < 0.0 ) THEN + WRITE(*,*) 'Warning! ', r_sq + r_sq = 0.0 + END IF + r_g = r_g + SQRT(r_sq) ! Accumulate root-mean-square radius of gyration + END DO ! End loop over polymers + + r_g = box * r_g / REAL(n) ! Average RMS Rg in sigma=1 units + + END FUNCTION rad_gyr + END PROGRAM qmc_pi_lj diff --git a/qmc_pi_lj_module.f90 b/qmc_pi_lj_module.f90 index a4e8ab9..9d33114 100644 --- a/qmc_pi_lj_module.f90 +++ b/qmc_pi_lj_module.f90 @@ -45,6 +45,7 @@ MODULE qmc_module ! Public derived type TYPE, PUBLIC :: potential_type ! A composite variable for interaction energies comprising REAL :: pot ! the cut (but not shifted) potential energy and + REAL :: vir ! the virial and LOGICAL :: ovr ! a flag indicating overlap (i.e. pot too high to use) CONTAINS PROCEDURE :: add_potential_type @@ -60,6 +61,7 @@ FUNCTION add_potential_type ( a, b ) RESULT (c) TYPE(potential_type) :: c ! Result is the sum of CLASS(potential_type), INTENT(in) :: a, b ! the two inputs c%pot = a%pot + b%pot + c%vir = a%vir + b%vir c%ovr = a%ovr .OR. b%ovr END FUNCTION add_potential_type @@ -68,6 +70,7 @@ FUNCTION subtract_potential_type ( a, b ) RESULT (c) TYPE(potential_type) :: c ! Result is the difference of CLASS(potential_type), INTENT(in) :: a, b ! the two inputs c%pot = a%pot - b%pot + c%vir = a%vir - b%vir c%ovr = a%ovr .OR. b%ovr ! This is meaningless, but inconsequential END FUNCTION subtract_potential_type @@ -114,11 +117,12 @@ END SUBROUTINE deallocate_arrays FUNCTION potential ( box, r_cut ) RESULT ( total ) IMPLICIT NONE - TYPE(potential_type) :: total ! Returns a composite of pot and ovr + TYPE(potential_type) :: total ! Returns a composite of pot, vir and ovr REAL, INTENT(in) :: box ! Simulation box length REAL, INTENT(in) :: r_cut ! Potential cutoff distance ! total%pot is the nonbonded cut (not shifted) classical potential energy for whole system + ! total%vir is the corresponding virial for whole system ! total%ovr is a flag indicating overlap (potential too high) to avoid overflow ! If this flag is .true., the values of total%pot etc should not be used ! Actual calculation is performed by function potential_1 @@ -135,7 +139,7 @@ FUNCTION potential ( box, r_cut ) RESULT ( total ) STOP 'Error in potential' END IF - total = potential_type ( pot=0.0, ovr=.FALSE. ) ! Initialize + total = potential_type ( pot=0.0, vir=0.0, ovr=.FALSE. ) ! Initialize DO k = 1, p ! Loop over ring polymers DO i = 1, n - 1 ! Loop over atoms within polymer @@ -158,7 +162,7 @@ END FUNCTION potential FUNCTION potential_1 ( rik, i, k, box, r_cut, j_range ) RESULT ( partial ) IMPLICIT NONE - TYPE(potential_type) :: partial ! Returns a composite of pot and ovr for given atom + TYPE(potential_type) :: partial ! Returns a composite of pot, vir and ovr for given atom REAL, DIMENSION(3), INTENT(in) :: rik ! Coordinates of atom of interest INTEGER, INTENT(in) :: i, k ! Index and polymer id of atom of interest REAL, INTENT(in) :: box ! Simulation box length @@ -166,6 +170,7 @@ FUNCTION potential_1 ( rik, i, k, box, r_cut, j_range ) RESULT ( partial ) INTEGER, OPTIONAL, INTENT(in) :: j_range ! Optional partner index range ! partial%pot is the nonbonded cut (not shifted) classical potential energy of atom rik with a set of other atoms + ! partial%vir is the corresponding virial of atom rik ! partial%ovr is a flag indicating overlap (potential too high) to avoid overflow ! If this is .true., the value of partial%pot should not be used ! The coordinates in rik are not necessarily identical with those in r(:,i,k) @@ -177,7 +182,7 @@ FUNCTION potential_1 ( rik, i, k, box, r_cut, j_range ) RESULT ( partial ) INTEGER :: j, j1, j2 REAL :: r_cut_box, r_cut_box_sq, box_sq - REAL :: sr2, sr6, r_ik_jk_sq + REAL :: sr2, sr6, sr12, r_ik_jk_sq REAL, DIMENSION(3) :: r_ik_jk REAL, PARAMETER :: sr2_ovr = 1.77 ! overlap threshold (pot > 100) TYPE(potential_type) :: pair @@ -212,7 +217,7 @@ FUNCTION potential_1 ( rik, i, k, box, r_cut, j_range ) RESULT ( partial ) r_cut_box_sq = r_cut_box**2 box_sq = box**2 - partial = potential_type ( pot=0.0, ovr=.FALSE. ) ! Initialize + partial = potential_type ( pot=0.0, vir=0.0, ovr=.FALSE. ) ! Initialize DO j = j1, j2 ! Loop over selected range of partners @@ -234,7 +239,9 @@ FUNCTION potential_1 ( rik, i, k, box, r_cut, j_range ) RESULT ( partial ) END IF sr6 = sr2**3 - pair%pot = sr6**2 - sr6 ! LJ potential (cut but not shifted) + sr12 = sr6**2 + pair%pot = sr12 - sr6 ! LJ potential (cut but not shifted) + pair%vir = pair%pot + sr12 ! LJ pair virial partial = partial + pair @@ -244,6 +251,7 @@ FUNCTION potential_1 ( rik, i, k, box, r_cut, j_range ) RESULT ( partial ) ! Include numerical factors partial%pot = partial%pot * 4.0 / REAL ( p ) ! Classical potentials are weaker by a factor p + partial%vir = partial%vir * 24.0 / REAL(3*p) ! 24*epsilon and divide virial by 3 & by p partial%ovr = .FALSE. ! No overlaps detected (redundant, but for clarity) END FUNCTION potential_1 @@ -257,7 +265,7 @@ FUNCTION spring ( box, k_spring ) RESULT ( total ) ! total is the quantum spring potential for whole system ! Actual calculation is performed by subroutine spring_1 - INTEGER :: i, k + INTEGER :: i, k, kp REAL :: partial_pot IF ( n /= SIZE(r,dim=2) ) THEN ! should never happen @@ -271,87 +279,45 @@ FUNCTION spring ( box, k_spring ) RESULT ( total ) total = 0.0 - DO k = 1, p ! Loop over ring polymers - DO i = 1, n ! Loop over atoms within polymer + DO k = 1, p ! Loop over atoms within polymer + + kp = k+1 + IF ( kp > p ) kp = 1 + + DO i = 1, n ! Loop over ring polymers - partial_pot = spring_1 ( r(:,i,k), i, k, box, k_spring, gt ) + partial_pot = spring_1 ( r(:,i,k), r(:,i,kp), box, k_spring ) total = total + partial_pot - END DO ! End loop over atoms within polymer - END DO ! End loop over ring polymers + END DO ! End loop over ring polymers + + END DO ! End loop over atoms within polymer END FUNCTION spring - FUNCTION spring_1 ( rik, i, k, box, k_spring, l_range ) RESULT ( partial_pot ) + FUNCTION spring_1 ( rik, ril, box, k_spring ) RESULT ( partial_pot ) IMPLICIT NONE REAL :: partial_pot ! Returns quantum potential for given atom REAL, DIMENSION(3), INTENT(in) :: rik ! Coordinates of atom of interest - INTEGER, INTENT(in) :: i, k ! Index and polymer id of atom of interest + REAL, DIMENSION(3), INTENT(in) :: ril ! Coordinates of other atom of interest REAL, INTENT(in) :: box ! Simulation box length REAL, INTENT(in) :: k_spring ! Spring potential constant - INTEGER, OPTIONAL, INTENT(in) :: l_range ! Optional partner index range - ! partial_pot is the quantum spring energy of atom (i,k) in rik for polymer k - ! with its neighbours (i,k-1) and (i,k+1), in the ring 1..p + ! partial_pot is the quantum spring energy of atom (i,k) in rik for polymer bead k + ! with its neighbour (i,l) in the same polymer ring ! The coordinates in rik are not necessarily identical with those in r(:,i,k) ! The partner atom index is always the same as i - ! The optional argument l_range restricts partner polymer ids to l=k-1, or l=k+1 ! It is assumed that r has been divided by box ! Results are in LJ units where sigma = 1, epsilon = 1 - INTEGER :: l, l1, l2 - REAL :: box_sq REAL :: r_ik_il_sq REAL, DIMENSION(3) :: r_ik_il - IF ( n /= SIZE(r,dim=2) ) THEN ! should never happen - WRITE ( unit=error_unit, fmt='(a,2i15)' ) 'Array bounds error for r, n=', n, SIZE(r,dim=2) - STOP 'Error in spring_1' - END IF - IF ( p /= SIZE(r,dim=3) ) THEN ! should never happen - WRITE ( unit=error_unit, fmt='(a,2i15)' ) 'Array bounds error for r, p=', p, SIZE(r,dim=3) - STOP 'Error in spring_1' - END IF - - IF ( PRESENT ( l_range ) ) THEN - SELECT CASE ( l_range ) - CASE ( lt ) ! Look at l = k-1 only - l1 = k-1 - l2 = k-1 - CASE ( gt ) ! Look at l = k+1 only - l1 = k+1 - l2 = k+1 - CASE default ! should never happen - WRITE ( unit = error_unit, fmt='(a,i10)') 'l_range error ', l_range - STOP 'Impossible error in spring_1' - END SELECT - ELSE ! Look at both l = k-1 and k+1 - l1 = k-1 - l2 = k+1 - END IF - - box_sq = box**2 - partial_pot = 0.0 ! Initialize - - DO l = l1, l2, 2 ! Loop over neighbours skipping l=k - - IF ( l < 1 ) THEN - r_ik_il(:) = rik(:) - r(:,i,p) ! link to end of ring polymer - ELSE IF ( l > p ) THEN - r_ik_il(:) = rik(:) - r(:,i,1) ! link to start of ring polymer - ELSE - r_ik_il(:) = rik(:) - r(:,i,l) - END IF - - r_ik_il(:) = r_ik_il(:) - ANINT ( r_ik_il(:) ) ! Periodic boundaries in box=1 units - r_ik_il_sq = SUM ( r_ik_il**2 ) * box_sq ! Squared distance in sigma=1 units - partial_pot = partial_pot + r_ik_il_sq ! Spring potential - - END DO ! End loop over possibilities skipping l=k - - ! Include numerical factors - partial_pot = partial_pot * 0.5 * k_spring + r_ik_il(:) = rik(:) - ril(:) + r_ik_il(:) = r_ik_il(:) - ANINT ( r_ik_il(:) ) ! Periodic boundaries in box=1 units + r_ik_il_sq = SUM ( r_ik_il**2 ) * box**2 ! Squared distance in sigma=1 units + partial_pot = 0.5 * k_spring * r_ik_il_sq ! Spring potential END FUNCTION spring_1