diff --git a/Examples/BK3.dat b/Examples/BK3.dat new file mode 100644 index 00000000..ad092b1a --- /dev/null +++ b/Examples/BK3.dat @@ -0,0 +1,10 @@ +1 +298.0 0.0 +O1 1.000000 1.00000 0.000 !oxygen +DM 1.162874 1.21043 0.000 !virtual site +G1 1.162874 1.21043 0.000 !drude particle attached to the virtual site +H2 1.975000 1.00000 0.000 !hydrogen 1 +G2 1.975000 1.00000 0.000 !drude particle attached to hydrogen 1 +H3 0.755550 1.94386 0.000 !hydrogen 2 +G3 0.755550 1.94386 0.000 !drude particle attached to hydogen 2 + diff --git a/Examples/BK3.mcf b/Examples/BK3.mcf new file mode 100644 index 00000000..ed8f21a6 --- /dev/null +++ b/Examples/BK3.mcf @@ -0,0 +1,73 @@ +!******************************************************************************** +! Cassandra - An open source atomistic Monte Carlo software package +! developed at the University of Notre Dame. +! http://cassandra.nd.edu +! Prof. Edward Maginn +! Copyright (2013) University of Notre Dame du Lac +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see . +!******************************************************************************** +!******************************************************************************* +! +! This is an input file for water BK3 +! +! Force field - BK3 +! J. Chem. Phys. 2013, 138, 204507 +!******************************************************************************* + +!********************************************************************************** +! +!********************************************************************************** + +! type element mass charge Exp-6 A (Kelvin) B (A^-1) C (Kelvin) rmin sigma alpha +# Atom_Info +7 +1 O1 O 16.00 0.0000 Born 38727700.7 3.56 399304.244 0.87 1.0 0.0 +2 DM D 0.000 0.0000 Born 0.000000 0.01 0.0 1.0 1.0 0.0 +3 G1 G 0.000 -1.168 Born 0.0 0.01 0.0 1.0 0.71 0.72 +4 H2 H 1.000 0.0000 Born 0.0 0.01 0.0 1.0 1.00 0.0 +5 G2 G 0.000 0.584 Born 0.0 0.01 0.0 1.0 0.40 0.36 +6 H3 H 1.000 0.0000 Born 0.0 0.01 0.0 1.0 1.00 0.0 +7 G3 G 0.000 0.584 Born 0.0 0.01 0.0 1.0 0.40 0.36 + + +# Bond_Info +6 +1 1 2 fixed 0.2661 +2 2 3 harmonic 131621.3 0.000 !0.5*spring constant for the drude particle (K = q*q/alpha), equilibrium bond length l=0 (shell rests on core) +3 1 4 fixed 0.9750 +4 1 6 fixed 0.9750 +5 4 5 harmonic 65810.64 0.000 +6 6 7 harmonic 65810.64 0.000 + +# Angle_Info +3 +1 4 1 6 fixed 104.52 +2 2 1 4 fixed 52.26 +3 2 1 6 fixed 52.26 + +# Dihedral_Info +0 + +# Fragment_Info +1 +1 7 1 2 3 4 5 6 7 + +# Improper_Info +0 + +# Fragment_Connectivity +0 + +END diff --git a/Examples/BK3_Cl.dat b/Examples/BK3_Cl.dat new file mode 100644 index 00000000..49495afb --- /dev/null +++ b/Examples/BK3_Cl.dat @@ -0,0 +1,4 @@ +1 +298.0 0.0 +Cl 1.00 2.00 0.000 !Chloride (core) +G1 1.00 2.00 0.000 !drude particle (shell) diff --git a/Examples/BK3_Na.dat b/Examples/BK3_Na.dat new file mode 100644 index 00000000..ff820d46 --- /dev/null +++ b/Examples/BK3_Na.dat @@ -0,0 +1,6 @@ +1 +298.0 0.0 +Na 1.50 1.20 0.000 !Na atom (core) +G1 1.50 1.20 0.000 !drude particle (shell) + + diff --git a/Examples/Cl_BK3.mcf b/Examples/Cl_BK3.mcf new file mode 100644 index 00000000..b7c90ac9 --- /dev/null +++ b/Examples/Cl_BK3.mcf @@ -0,0 +1,59 @@ +!******************************************************************************** +! Cassandra - An open source atomistic Monte Carlo software package +! developed at the University of Notre Dame. +! http://cassandra.nd.edu +! Prof. Edward Maginn +! Copyright (2013) University of Notre Dame du Lac +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see . +!******************************************************************************** +!******************************************************************************* +! +! This is an input file for BK3 Cl. +! + +!******************************************************************************* + +!********************************************************************************** +!Molecular connectivity file for water.pdb +!********************************************************************************** + + +# Atom_Info +2 +1 Cl Cl 35.45 -11.0 Born 113056021.0 3.10 962178.90 0.788 0.96 0.0 +2 G1 G 0.000 10.0 Born 0.0 0.01 0.0 1.000 0.96 3.5 + + +# Bond_Info +1 +1 1 2 harmonic 1984745.7 0.0 + +# Angle_Info +0 + +# Dihedral_Info +0 + +# Fragment_Info +1 +1 2 1 2 + +# Improper_Info +0 + +# Fragment_Connectivity +0 + +END diff --git a/Examples/H2O_NaCl_GEMC.inp b/Examples/H2O_NaCl_GEMC.inp new file mode 100644 index 00000000..9d3c72f3 --- /dev/null +++ b/Examples/H2O_NaCl_GEMC.inp @@ -0,0 +1,151 @@ +!****************************************************************************** +! +! This is an input file for SPC/E water. +! +! Created by Eliseo Marin 12/13/13 +!******************************************************************************* + + +# Run_Name +BK3_salt_H2O_GEMC +!-------------- + +# Sim_Type +GEMC +!--------------- + +# Nbr_Species +3 +!---------------- + +# VDW_Style +Born cut_tail 12.0 !Born (bukcingham potential) cut-off distance 12.0 A +Born cut_tail 19.0 +!---------------- + +# Rcutoff_Low +0.1 +!---------------- + +# Mixing_Rule +kong ! Kong mixing rule for buckingham potential +!---------------- + +# Charge_Style +coul Gaussian_charge 12.00 0.0001 !Ewald sum for Gaussian charge +coul Gaussian_charge 19.00 0.0001 +!---------------- + +# Intra_Scaling +0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 +!------------------ + +# Molecule_Files +BK3.mcf 256 +Na_BK3.mcf 2 +Cl_BK3.mcf 2 +!--------------------- + +# Pair_Energy +FALSE !Pair_Energy must be set to FALSE + +# Box_Info +2 +CUBIC +30.0 + +CUBIC +40.0 + +# Temperature_Info +473.15 473.15 +!------------------ + +# Seed_Info +2082012 82012 +!------------------- + +# Move_Probability_Info + +# Prob_Translation +0.25 +0.002 0.002 0.002 !adjustable parameter for multiparticle translation move +0.002 0.002 0.002 + +# Prob_Rotation +0.25 +0.001 0.001 0.001 !adjustable parameter for multiparticle rotation move +0.001 0.001 0.001 + +# Prob_Volume +0.1 +300 +300 + +# Prob_Swap +0.4 +insertion method +reservior +insertion method +none +insertion method +none +# Prob_Species_Swap +1.0 1.0 1.0 + +# Done_Probability_Info + +!---------------------- +# Bond_Prob_Cutoff +1.0E-10 + +# CBMC_Info +kappa_ins 10 +kappa_rot 0 +kappa_dih 10 +rcut_cbmc 9.0 9.0 + +# Start_Type +make_config +250 6 +2 0 +2 0 + +# Run_Type +Equilibration 100 100 +!------------------ + +# Average_Info ( 0 == yes, 1 = No) +1 + +# Frequency_Info +freq_type none +Nthermofreq 100 +Ncoordfreq 1000 +MCsteps 1000 +# Done_Frequency_Info +!-------------------- + +# Property_Info 1 +Energy_Total +Volume +Nmols + +# Property_Info 2 +Energy_Total +Volume +Nmols +Pressure1 !Calculation of pressure by performing a (positive) trial volume change +Pressure2 !Calculation of pressure by performing a (negative) trial volume change (without evaluation of virial pressure) + +# Fragment_Files +BK3.dat 1 +BK3_Na.dat 2 +BK3_Cl.dat 3 + +END diff --git a/Examples/H2O_NaCl_NPT.inp b/Examples/H2O_NaCl_NPT.inp new file mode 100644 index 00000000..6a1ec0b4 --- /dev/null +++ b/Examples/H2O_NaCl_NPT.inp @@ -0,0 +1,128 @@ +!****************************************************************************** +! +! This is an input file for SPC/E water. +! +! Created by Eliseo Marin 12/13/13 +!******************************************************************************* + + +# Run_Name +BK3_salt_H2O_NPT +!-------------- + +# Sim_Type +NPT_MC +!--------------- + +# Nbr_Species +3 +!---------------- + +# VDW_Style +Born cut_tail 10.0 !Born (bukcingham potential) cut-off distance 12.0 A +!---------------- + +# Rcutoff_Low +0.1 +!---------------- + +# Mixing_Rule +kong ! Kong mixing rule for buckingham potential +!---------------- + +# Charge_Style +coul Gaussian_charge 10.00 0.0001 !Ewald sum for Gaussian charge +!---------------- + +# Intra_Scaling +0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 +!------------------ + +# Molecule_Files +BK3.mcf 256 +Na_BK3.mcf 2 +Cl_BK3.mcf 2 +!--------------------- + +# Pair_Energy +FALSE !Pair_Energy must be set to FALSE + +# Box_Info +1 +CUBIC +30.0 + +# Temperature_Info +298.15 + +# Pressure_Info +1.0 +!------------------ + +# Seed_Info +2082012 82012 +!------------------- + +# Move_Probability_Info + +# Prob_Translation +0.3 +0.002 0.002 0.002 !adjustable parameter for multiparticle translation move + +# Prob_Rotation +0.3 +0.001 0.001 0.001 !adjustable parameter for multiparticle rotation move + +# Prob_Volume +0.4 +300 + +# Done_Probability_Info + +!---------------------- +# Bond_Prob_Cutoff +1.0E-10 + +# CBMC_Info +kappa_ins 10 +kappa_rot 0 +kappa_dih 10 +rcut_cbmc 9.0 9.0 + +# Start_Type +make_config +256 +2 +2 + +# Run_Type +Equilibration 100 100 +!------------------ + +# Average_Info ( 0 == yes, 1 = No) +1 + +# Frequency_Info +freq_type none +Nthermofreq 100 +Ncoordfreq 100 +MCsteps 500 +# Done_Frequency_Info +!-------------------- + +# Property_Info 1 +Energy_Total +Volume +Density + +# Fragment_Files +BK3.dat 1 +BK3_Na.dat 2 +BK3_Cl.dat 3 + +END diff --git a/Examples/Na_BK3.mcf b/Examples/Na_BK3.mcf new file mode 100644 index 00000000..f7aaafdc --- /dev/null +++ b/Examples/Na_BK3.mcf @@ -0,0 +1,52 @@ +!******************************************************************************** +! Cassandra - An open source atomistic Monte Carlo software package +! developed at the University of Notre Dame. +! http://cassandra.nd.edu +! Prof. Edward Maginn +! Copyright (2013) University of Notre Dame du Lac +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see . +!******************************************************************************** +!******************************************************************************* +! +!******************************************************************************* + + +# Atom_Info +2 +1 Na Na 23.0 11.0 Born 1.02231508E+015 15.0 66149.79 0.03 0.65 0.0 +2 G1 G 0.00 -10.0 Born 0.0 0.01 0.0 1.00 0.65 0.157 + + +# Bond_Info +1 +1 1 2 harmonic 44245923.0 0.0 + +# Angle_Info +0 + +# Dihedral_Info +0 + +# Fragment_Info +1 +1 2 1 2 + +# Improper_Info +0 + +# Fragment_Connectivity +0 + +END diff --git a/Examples/polarizable.tex b/Examples/polarizable.tex new file mode 100644 index 00000000..509f91a8 --- /dev/null +++ b/Examples/polarizable.tex @@ -0,0 +1,163 @@ +\documentclass[12pt]{article} +\usepackage[latin1]{inputenc} +\usepackage[reqno]{amsmath} +\usepackage{amssymb,latexsym} +\usepackage[pdftex]{color,graphicx} +\usepackage{multirow} +\usepackage{multicol} +\usepackage{caption} +\usepackage{float} +%\usepackage{cite} +\usepackage{gensymb} +\usepackage{color} \newcommand{\cor}{\color{red}} +\usepackage[sort, numbers, super, compress]{natbib} +\usepackage{indentfirst} +\setlength{\parindent}{24pt} +\renewcommand{\thesection}{\Roman{section}.} +\renewcommand{\thesubsection}{\Alph{subsection}.} + +\title{\textbf{Simulation of Drude oscillator based polarizable force field with Cassandra}} +\date{} +\begin{document} +\maketitle +\section{Algorithm} +\subsection{Multi-particle move} +Force/torque biased multi-particle translation/rotation moves were implemented (mpm\_translate\_rotate.f90). All particles are translated or rotated simultaneously in one step. For translation move, the translation vector for particle $k$ is given as, +\begin{equation} +\boldsymbol{t}_{k} = \beta A \boldsymbol{F}_{k}^{old} + \boldsymbol{\gamma}_{k} +\end{equation} +\noindent{$\boldsymbol{F}_{k}^{old}$ is the force on particle $k$, $\boldsymbol{\gamma}_{k}$ is a random number chosen from a Gaussian distribution with 0 mean and variance $2 A$, and $A$ is adjusted during the simulation to reach a certain acceptance ratio (0.33).} The acceptance criteria of the multi-particle translation move is given by,\\ +\begin{equation} + s^{old->new} = \sum_{k=1}^{N} -(\beta A \boldsymbol{F}_{k}^{old} - \boldsymbol{\gamma}_{k} )^2/(4A) +\end{equation} +\begin{equation} + s^{new->old} = \sum_{k=1}^{N} -(\beta A \boldsymbol{F}_{k}^{new} + \boldsymbol{\gamma}_{k} )^2/(4A) +\end{equation} +\begin{equation} + \mathrm{Prob} = min \big( 1,\mathrm{exp}(-\beta \Delta U + s^{new->old} - s^{old->new} \big) +\end{equation} +The details for the multi-particle translation/rotation move can be found in Moucka et al., {\it Mol. Sim.} 2013, 39, 1125-1134. +\subsection{Ewald summation for Gaussian charge} + +A Gaussian charge is represented by a spherical charge distribution, +\begin{equation} + \rho_{i}(r) = \frac{q_i}{(2 \pi \sigma_i^2)^{3/2}} \mathrm{exp}\Big(\frac{-|\boldsymbol{r}-\boldsymbol{r_i}|}{2\sigma_i^2}\Big) +\end{equation} +\noindent{where $\sigma_i$ is the with of the charge distribution. The interaction between two Gaussian charges is given by,} +\begin{equation} + U_{ij} = \frac{q_i q_j}{r_{ij}} \mathrm{erf}(\alpha_{ij} r_{ij}) +\end{equation} +\noindent{where $\alpha_{ij} =1/ \sqrt{2 (\sigma_i^2 + \sigma_j^2)} $.} The direct-space part of the Ewald summation is given as, +\begin{equation} + U^{dir} = \sum_{i=1}^{N} \sum_{j=i+1}^{N} \frac{q_i q_j}{r_{ij}} \Big[ \mathrm{erf}(\alpha_{ij}r_{ij}) - \mathrm{erf}(k r_{ij}) \Big] +\end{equation} +\noindent{$k$ is the convergence parameter for the Ewald summation, and the reciprocal part of the summation is the same as that for point charge. The details of the calculation can be found at Kiss et al., {\it J. Chem. Theory. Comput.} 2014, 10, 5513-5519.} \\ + +\subsection{Energy minimization for Drude particles} +The forces acting on Drude particles are calculated at every MC step, and the positions of Drude particles are calculated by the following iteration, +\begin{equation} + \boldsymbol{r}_{kD}(n) = \boldsymbol{r}_{kD}(n-1) + \boldsymbol{F}_{kD}/K +\end{equation} +\noindent{where $K$ is the spring constant ($K = q^2/\alpha$, $\alpha$ is the polarizability). The iteration is terminated if +\begin{equation} + \max_{i=1...N}|\boldsymbol{r}_{kD}(n) - \boldsymbol{r}_{kD}(n-1)| < 10^{-4} \mathrm{nm} +\end{equation} + + +\section{Simulation Input} + +A set of recently proposed polarizable H$_2$O, Na$^+$ and Cl$^-$ force field models (BK3 models) is used here as an example. The vdW interactions of the BK3 force fields are represented by the Buckingham Exp-6 potential. The details of the force fields can be found in {\it J. Chem. Phys.} 2013, 138, 204507 and 2014, 141, 114501. + +\subsection{MCF File} +\noindent{\# Atom\_Info} \\ +Integer(1) \\ +Integer(2) Character(3) Character(4) Real(5) Real(6) Character(7) Real(8) Real(9) Real(10) Real(11) Real(12) Real(13) \\ + +\noindent{Integer(1): total number of atoms in the species.} \\ + +\noindent{Integer(2): index of the atom.} \\ + +\noindent{Character(3): type of the atom.} \\ + +\noindent{Character(4): element name of the atom. `G' indicates that the atom is a Drude particle (shell). The Drude particle (shell) must follow its connecting atom (core) in MCF file.} \\ + +\noindent{Real(5): mass of the atom.} \\ + +\noindent{Real(6): charge (e).} \\ + +\noindent{Character(7): function form for vdW interaction. `Born' is the Born (Buckingham Exp-6) potential. +$U = A\mathrm{exp}(-B r) - C^6/r^6$.} \\ + +\noindent{Real(8): A parameter of the Exp-6 potential, in Kelvin.} \\ + +\noindent{Real(9): B parameter of the Exp-6 potential, in \AA$^{-1}$.} \\ + +\noindent{Real(10): C parameter of the Exp-6 potential, in Kelvin.} \\ + +\noindent{Real(10): minimum cut off distance for the Exp-6 potential, in \AA. } \\ + +\noindent{Real(11): width of the Gaussian charge ($\sigma$), in \AA.} \\ + +\noindent{Real(12): polarizablility of the atom ($\alpha$), in \AA$^3$.} \\ + +\noindent{\# Bond\_Info} \\ +Integer(1) \\ +Integer(2) Integer(3) Integer(4) Character(5) Real(6) Real(7) \\ + +\noindent{Integer(1): total number of bonds in the species.} \\ + +\noindent{Integer(2): index of the bond.} \\ + +\noindent{Integer(3), Integer(4): IDs of the atoms in the bond.} \\ + +\noindent{Character(5): type of the bond. `fixed' for fixed bond and `harmonic' for harmonic bond. For a core-shell unit, the bond type must be `harmonic' and the spring constant ($\frac{1}{2}q^2/\alpha$) must be provided.} \\ + +\noindent{Real(6): bond length in $\AA$ if the bond type is `fixed' or spring constant ($\frac{1}{2}K)$ if the bond type is 'harmonic'.} \\ + +\noindent{Real(7): bond length for harmonic bond in $\AA$. For a core-shell unit, this value is 0.} + +\subsection{Input File} + +\noindent{\# VDW\_Style} \\ +Character(1) Character(2) Real(3) \\ + +\noindent{Character(1): type of vdW interaction. `LJ' for Lennard-Jones potential and `Born' for Born (Buckingham exp-6) potential.} \\ + +\noindent{Character(2): cut-off scheme for vdW interaction. For Buckingham potential, only `cut\_tail' is permitted.} \\ + +\noindent{Real(3): cut off distance in \AA.} \\ + +\noindent{\# Mixing\_Rule} \\ +Character(1) \\ + +\noindent{Character(1): `LB', `Geometric' and `Kong' are permitted. `Kong' is the mixing rule for Buckingham potential.} \\ + +\noindent{\# Pair\_Energy} \\ +Character(1) \\ + +\noindent{Character(1): For simulation of polarizable model, use `FALSE'.} \\ + + +\noindent{\# Prob\_Translation} \\ +Real(1) \\ +Real(i,j) \\ + +\noindent{Real(1): Probability for the move.} \\ + +\noindent{Real(i,j): Adjustable parameter for the multi-particle move or maximum displacement for the single particle move.} \\ + + +\noindent{\# Property\_Info} \\ + +`Pressure1': calculate pressure by the energy fluctuation due to a positive (volume of the box increases) trial volume change. \\ + +`Pressure2': calculate pressure by the energy fluctuation due to a negative (volume of the box decreases) trial volume change. \\ + +The pressure of a system is estimated by ({\it J. Chem. Phys., 1996, 105, 8469}), +\begin{equation} +P = \frac{k_B T}{\Delta V}\mathrm{ln}\big<\frac{V'}{V}\mathrm{exp}(-\beta \Delta U)\big> +\end{equation} + +\noindent{$\Delta V$ is set to 25 $\AA^3$, $V'=V\pm\Delta V$, and the instantaneous value of $\frac{V'}{V}\mathrm{exp}(-\beta \Delta U)$ is output. } + +\end{document} diff --git a/Src/Makefile b/Src/Makefile index 4c6c325d..3204cb86 100644 --- a/Src/Makefile +++ b/Src/Makefile @@ -1,3 +1,4 @@ + #!******************************************************************************** #! Cassandra - An open source atomistic Monte Carlo software package #! developed at the University of Notre Dame. @@ -61,8 +62,8 @@ FC = ifort LIBS = INCS = # -#F90FLAGS = -c -ffree-line-length-none -F90FLAGS = -c -g -check all -traceback #-warn unused #-openmp +F90FLAGS = -c -ffree-line-length-none +#F90FLAGS = -c -g -check all -traceback #-warn unused #-openmp #F90FLAGS = -c -openmp #F90FLAGS = -p -c -O3 #F90FLAGS = -c -O3 -openmp @@ -128,6 +129,10 @@ OBJS = main.o \ atom_displacement.o \ chempot.o \ pair_nrg_variables.o \ + shell_gemc_nvt_volume.o\ + shell_gemc_particle_transfer.o\ + shell_relax.o\ + mpm_translate_rotate.o #################################################################### diff --git a/Src/Makefile.intel.openMP b/Src/Makefile.intel.openMP index 61df9ea0..e857274f 100644 --- a/Src/Makefile.intel.openMP +++ b/Src/Makefile.intel.openMP @@ -1,3 +1,4 @@ + #!******************************************************************************** #! Cassandra - An open source atomistic Monte Carlo software package #! developed at the University of Notre Dame. @@ -57,7 +58,7 @@ FC = ifort LIBS = INCS = # -F90FLAGS = -c -O3 -openmp +F90FLAGS = -O -c -openmp LINKFLAGS = -openmp @@ -119,6 +120,12 @@ OBJS = main.o \ atom_displacement.o \ chempot.o \ pair_nrg_variables.o \ + shell_gemc_nvt_volume.o\ + shell_gemc_particle_transfer.o\ + shell_relax.o\ + mpm_translate_rotate.o + + #################################################################### diff --git a/Src/accept_or_reject.f90 b/Src/accept_or_reject.f90 index 6dd2e5a7..e9b355d4 100644 --- a/Src/accept_or_reject.f90 +++ b/Src/accept_or_reject.f90 @@ -1,3 +1,4 @@ + !******************************************************************************** ! Cassandra - An open source atomistic Monte Carlo software package ! developed at the University of Notre Dame. diff --git a/Src/create_nonbond_table.f90 b/Src/create_nonbond_table.f90 index 27fed2dd..1784a789 100644 --- a/Src/create_nonbond_table.f90 +++ b/Src/create_nonbond_table.f90 @@ -1,55 +1,62 @@ !******************************************************************************** -! CASSANDRA - Computational Atomistic Simulation Software at Notre Dame -! for Research in Academia. -! http://molsim.wiki.zoho.com/ -! Copyright (2007) University of Notre Dame. -! Authors: Ed Maginn (ed@nd.edu); Jindal Shah (jshah@nd.edu) -!******************************************************************************** - -!******************************************************************************** - SUBROUTINE Create_Nonbond_Table -!******************************************************************************** - ! This routine determines the number and identity of all unique atom types - ! in the system. It then takes the associated vdw paramters of each atom - ! and creates an interaction table using the specified mixing rule. - ! Currently only supports LJ vdw type and LB or geometric mixing rules. - - ! The final product of this code are a bunch of matrices - ! vdw_param1_table(itype,jtype), vdw_param2_table(itype,jtype), etc. - ! where for the LJ potential, vdw_param1_table(itype,jtype) is the epsilon - ! to use between atoms of itype and jtype, and vdw_param2_table(itype,jtype) - ! is the sigma to use between these two atom type. - - ! Example: let's say you need to compute the LJ energy betweem atom 4 of species 2 - ! and atom 6 of species 1. - ! You first identify the unique atom type of these two atoms: - ! itype = nonbond_list(4,2)%atom_type_number - ! jtype = nonbond_list(6,1)%atom_type_number - - ! Then get the parameters: - ! epsilon(itype, jtype) = vdw_param1_table(itype,jtype) - ! sigma(itype,jtype) = vdw_param2_table(itype,jtype) - ! Note that these matrices are symmetric, so - ! vdw_param1_table(itype,jtype) = vdw_param1_table(jtype,itype) - - ! Written: Sat Oct 6 08:58:06 MDT 2007 - ! Author: E. Maginn - - ! *** CALLS *** - ! Clean_Abort - ! - ! *** CALLED BY *** - ! NVTMC_Control - ! - ! Revision history: - - ! 04/08/09 (TR) : Removed any character evaluations. - ! - ! 08/26/11 (JS) : Steele potential constants calculations - ! - ! 03/27/14 Eliseo Rimoldi : custom (manually input) nonbond potentials added - -!******************************************************************************** +! Cassandra - An open source atomistic Monte Carlo software package +! developed at the University of Notre Dame. +! http://cassandra.nd.edu +! Prof. Edward Maginn +! Copyright (2013) University of Notre Dame du Lac +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see . +!******************************************************************************* + + +SUBROUTINE Create_Nonbond_Table + !******************************************************************************** + ! This routine determines the number and identity of all unique atom types + ! in the system. It then takes the associated vdw paramters of each atom + ! and creates an interaction table using the specified mixing rule. + ! Currently only supports LJ vdw type and LB or geometric mixing rules. + + ! The final product of this code are a bunch of matrices + ! vdw_param1_table(itype,jtype), vdw_param2_table(itype,jtype), etc. + ! where for the LJ potential, vdw_param1_table(itype,jtype) is the epsilon + ! to use between atoms of itype and jtype, and vdw_param2_table(itype,jtype) + ! is the sigma to use between these two atom type. + + ! Example: let's say you need to compute the LJ energy betweem atom 4 of species 2 + ! and atom 6 of species 1. + ! You first identify the unique atom type of these two atoms: + ! itype = nonbond_list(4,2)%atom_type_number + ! jtype = nonbond_list(6,1)%atom_type_number + + ! Then get the parameters: + ! epsilon(itype, jtype) = vdw_param1_table(itype,jtype) + ! sigma(itype,jtype) = vdw_param2_table(itype,jtype) + ! Note that these matrices are symmetric, so + ! vdw_param1_table(itype,jtype) = vdw_param1_table(jtype,itype) + ! + ! Called by + ! + ! gcmc_control + ! gemc_control + ! nptmc_control + ! nvtmc_control + ! nvt_mc_fragment_control + ! + ! Revision history + ! + ! 12/10/13 : Beta Release + !******************************************************************************** USE Run_Variables USE Type_Definitions USE IO_Utilities @@ -65,17 +72,13 @@ SUBROUTINE Create_Nonbond_Table ! Steele potential REAL(DP) :: sigma_ss, eps_ss, rho_s, delta_s, eps_sf - !custom mixing rules - INTEGER :: ierr,line_nbr,nbr_entries, is_1, is_2, ia_1, ia_2, itype_custom, jtype_custom - CHARACTER(120) :: line_string, line_array(20) - !******************************************************************************** ALLOCATE(temp_atomtypes(1000), Stat=AllocateStatus) IF (AllocateStatus .NE. 0) THEN err_msg = '' err_msg(1) = ' ERROR: Not enough memory for temp_atomtypes ' - CALL Clean_Abort(err_msg,'create_nonbond_table') + CALL Clean_Abort(err_msg,'ceate_nonbond_table') END IF ! Initialize the number of atom types and the temp_atomtypes @@ -169,6 +172,7 @@ SUBROUTINE Create_Nonbond_Table ALLOCATE(vdw_param3_table(nbr_atomtypes,nbr_atomtypes), Stat=AllocateStatus) ALLOCATE(vdw_param4_table(nbr_atomtypes,nbr_atomtypes), Stat=AllocateStatus) ALLOCATE(vdw_param5_table(nbr_atomtypes,nbr_atomtypes), Stat=AllocateStatus) + ! ALLOCATE(vdw_param6_table(nbr_atomtypes,nbr_atomtypes), Stat=AllocateStatus) ! Allocate memory for total number bead types in each box ALLOCATE(nint_beads(nbr_atomtypes,nbr_boxes)) @@ -197,168 +201,189 @@ SUBROUTINE Create_Nonbond_Table DO jtype = 1, nbr_atomtypes + ! Flags that are tripped when a particular atomtype parameter set is located + iset = 0 + jset = 0 - IF (mix_rule /= 'custom') THEN + DO is = 1, nspecies - ! Flags that are tripped when a particular atomtype parameter set is located - iset = 0 - jset = 0 + DO ia = 1, natoms(is) - DO is = 1, nspecies + IF (iset == 0) THEN + ! Search for atomtype i because the parameters have not been found - DO ia = 1, natoms(is) + IF (nonbond_list(ia,is)%atom_type_number == itype) THEN + ! This type of atom has been found. Grab its vdw parameters - IF (iset == 0) THEN - ! Search for atomtype i because the parameters have not been found + DO k=1, nbr_vdw_params(is) - IF (nonbond_list(ia,is)%atom_type_number == itype) THEN - ! This type of atom has been found. Grab its vdw parameters - - DO k=1, nbr_vdw_params(is) - - temp_param_i(k) = nonbond_list(ia,is)%vdw_param(k) + temp_param_i(k) = nonbond_list(ia,is)%vdw_param(k) - ENDDO + ENDDO - ! This atomtype parameters are now determined. Trip flag - iset = 1 + ! This atomtype parameters are now determined. Trip flag + iset = 1 - atom_type_list(itype) = nonbond_list(ia,is)%atom_name - - ENDIF + atom_type_list(itype) = nonbond_list(ia,is)%atom_name - ENDIF + ENDIF + + ENDIF - IF(jset == 0) THEN + IF(jset == 0) THEN - ! Now search for parameters of type j + ! Now search for parameters of type j - IF (nonbond_list(ia,is)%atom_type_number == jtype) THEN + IF (nonbond_list(ia,is)%atom_type_number == jtype) THEN - DO k=1, nbr_vdw_params(is) + DO k=1, nbr_vdw_params(is) - temp_param_j(k) = nonbond_list(ia,is)%vdw_param(k) + temp_param_j(k) = nonbond_list(ia,is)%vdw_param(k) - ENDDO + ENDDO - ! This atom type has been located - jset = 1 + ! This atom type has been located + jset = 1 - atom_type_list(jtype) = nonbond_list(ia,is)%atom_name + atom_type_list(jtype) = nonbond_list(ia,is)%atom_name - ENDIF + ENDIF - ENDIF + ENDIF - ENDDO + ENDDO - ENDDO + ENDDO - ! atom types i and j have been found and vdw parameters loaded into temporary arrays. - ! Apply mixing rules and load vdw_table + ! atom types i and j have been found and vdw parameters loaded into temporary arrays. + ! Apply mixing rules and load vdw_table - ! Adapt this to other potential types and mixing rules. Custom mixing rules can - ! be created manually and then this routine should be bypassed. + ! Adapt this to other potential types and mixing rules. Custom mixing rules can + ! be created manually and then this routine should be bypassed. - IF (int_vdw_style(1) == vdw_lj) THEN - ! There are two vdw parameters + IF (int_vdw_style(1) == vdw_lj) THEN + ! There are two vdw parameters - ! Set LJ epsilon - IF ( (temp_param_i(1) <= tiny_number) .OR. (temp_param_j(1) <= tiny_number) ) THEN - vdw_param1_table(itype,jtype) = 0.0_DP - ! for parameters with zero, avoid overflow and set to zero + ! Set LJ epsilon + IF ( (temp_param_i(1) <= tiny_number) .OR. (temp_param_j(1) <= tiny_number) ) THEN + vdw_param1_table(itype,jtype) = 0.0_DP + ! for parameters with zero, avoid overflow and set to zero - ELSE + ELSE - ! Use specified mixing rule + ! Use specified mixing rule - ! LB mixing rule: epsij = (epsi * epsj)^(1/2); sigmaij = 1/2 (sigmai + sigmaj) - IF (mix_rule == 'LB') & - vdw_param1_table(itype,jtype) = dsqrt(temp_param_i(1)*temp_param_j(1)) + ! LB mixing rule: epsij = (epsi * epsj)^(1/2); sigmaij = 1/2 (sigmai + sigmaj) + IF (mix_rule == 'LB') & + vdw_param1_table(itype,jtype) = dsqrt(temp_param_i(1)*temp_param_j(1)) - ! geometric mixing rule: epsij = (epsi * epsj)^(1/2); sigmaij = (sigmai * sigmaj)^(1/2) - IF (mix_rule == 'geometric') & - vdw_param1_table(itype,jtype) = dsqrt(temp_param_i(1)*temp_param_j(1)) + ! geometric mixing rule: epsij = (epsi * epsj)^(1/2); sigmaij = (sigmai * sigmaj)^(1/2) + IF (mix_rule == 'geometric') & + vdw_param1_table(itype,jtype) = dsqrt(temp_param_i(1)*temp_param_j(1)) - ENDIF + ENDIF + + ! Set LJ sigma + IF ( (temp_param_i(2) <= 1.0E-05) .OR. (temp_param_j(2) <= 1.0E-05) ) THEN + vdw_param2_table(itype,jtype) = 0.0_DP - ! Set LJ sigma - IF ( (temp_param_i(2) <= 1.0E-05) .OR. (temp_param_j(2) <= 1.0E-05) ) THEN - vdw_param2_table(itype,jtype) = 0.0_DP + ELSE - ELSE + IF (mix_rule == 'LB') & + vdw_param2_table(itype,jtype) = (temp_param_i(2) + temp_param_j(2)) * 0.5 + IF (mix_rule == 'geometric') & + vdw_param2_table(itype,jtype) = dsqrt(temp_param_i(2) * temp_param_j(2)) - IF (mix_rule == 'LB') & - vdw_param2_table(itype,jtype) = (temp_param_i(2) + temp_param_j(2)) * 0.5 - IF (mix_rule == 'geometric') & - vdw_param2_table(itype,jtype) = dsqrt(temp_param_i(2) * temp_param_j(2)) + ENDIF - ENDIF + ENDIF + + IF (int_vdw_style(1) == born_cut_tail) THEN + + IF ( (abs(temp_param_i(1)) <= tiny_number) .OR. (abs(temp_param_j(1)) <= tiny_number) ) THEN + vdw_param1_table(itype,jtype) = 0.0_DP + ! for parameters with zero, avoid overflow and set to zero + ELSE + ! Use specified mixing rule + ! LB mixing rule: epsij = (epsi * epsj)^(1/2); sigmaij = 1/2 (sigmai + sigmaj) + IF (mix_rule == 'LB') & + vdw_param1_table(itype,jtype) = dsqrt(temp_param_i(1)*temp_param_j(1)) + ! geometric mixing rule: epsij = (epsi * epsj)^(1/2); sigmaij = (sigmai * sigmaj)^(1/2) + IF (mix_rule == 'geometric') & + vdw_param1_table(itype,jtype) = dsqrt(temp_param_i(1)*temp_param_j(1)) + + IF (mix_rule == 'kong') THEN + vdw_param1_table(itype,jtype) = temp_param_i(1)*(temp_param_i(1)*temp_param_i(2)/temp_param_j(1)/temp_param_j(2))**(-temp_param_i(2)/(temp_param_i(2) + temp_param_j(2))) + & + & temp_param_j(1)*(temp_param_j(1)*temp_param_j(2)/temp_param_i(1)/temp_param_i(2))**(-temp_param_j(2)/(temp_param_i(2) + temp_param_j(2))) + vdw_param1_table(itype,jtype) = vdw_param1_table(itype,jtype) * 0.5 + END IF - ENDIF + ENDIF + + IF ( (abs(temp_param_i(2)) <= tiny_number) .OR. (abs(temp_param_j(2)) <= tiny_number) ) THEN + vdw_param2_table(itype,jtype) = 0.0_DP + ELSE + IF (mix_rule == 'LB') & + vdw_param2_table(itype,jtype) = (temp_param_i(2)+temp_param_j(2))/2.0 + IF (mix_rule == 'geometric') & + vdw_param2_table(itype,jtype) = dsqrt(temp_param_i(2)*temp_param_j(2)) + IF (mix_rule == 'kong') & + vdw_param2_table(itype,jtype) = 2.0*(temp_param_i(2)*temp_param_j(2)) / (temp_param_i(2)+temp_param_j(2)) + ENDIF + + IF ( (abs(temp_param_i(3)) <= tiny_number) .OR. (abs(temp_param_j(3)) <= tiny_number) ) THEN + vdw_param3_table(itype,jtype) = 0.0_DP + ELSE + IF (mix_rule == 'LB') & + vdw_param3_table(itype,jtype) = dsqrt(temp_param_i(3)*temp_param_j(3)) + IF (mix_rule == 'geometric') & + vdw_param3_table(itype,jtype) = dsqrt(temp_param_i(3)*temp_param_j(3)) + IF (mix_rule == 'kong') & + vdw_param3_table(itype,jtype) = dsqrt(temp_param_i(3)*temp_param_j(3)) + + ENDIF + + IF ( (abs(temp_param_i(4)) <= tiny_number) .OR. (abs(temp_param_j(4)) <= tiny_number) ) THEN + vdw_param4_table(itype,jtype) = 0.0_DP + ELSE + IF (mix_rule == 'LB') & + vdw_param4_table(itype,jtype) = (temp_param_i(4)+temp_param_j(4))/2.0 + IF (mix_rule == 'geometric') & + vdw_param4_table(itype,jtype) = dsqrt(temp_param_i(4)*temp_param_j(4)) + IF (mix_rule == 'kong') & + vdw_param4_table(itype,jtype) = 2.0*(temp_param_i(4)*temp_param_j(4))/(temp_param_i(4)+temp_param_j(4)) + ENDIF + + ENDIF - IF (int_vdw_style(1) == vdw_lj) THEN - ! Report parameters to logfile. Format is specific to vdw type. Add others here if - ! other than LJ potential is used. - - WRITE(logunit,'(A6,5x,A6,2x,T20,f10.4,T50,f10.4)') & - atom_type_list(itype), atom_type_list(jtype), & - vdw_param1_table(itype,jtype), vdw_param2_table(itype,jtype) - - ENDIF - - - ELSE !custom mixing rule if - - REWIND(inputunit) - - ierr = 0 - line_nbr = 0 - - DO - line_nbr = line_nbr + 1 - - CALL Read_String(inputunit,line_string,ierr) - - - IF (ierr .NE. 0) THEN - err_msg = "" - err_msg(1) = "Error reading mixing rules1." - CALL Clean_Abort(err_msg,'Get_Mixing_Rules') - END IF - - IF (line_string(1:13) == '# Mixing_Rule') THEN - READ(inputunit,*) - ! Assign the first entry on the line to the mixing rule - DO is_1 = 1, nspecies - DO ia_1 = 1, natoms(is_1) - itype_custom = nonbond_list(ia_1,is_1)%atom_type_number - DO is_2 = 1, nspecies - DO ia_2 = 1, natoms(is_2) - CALL Parse_String(inputunit,line_nbr,1,nbr_entries,line_array,ierr) - jtype_custom = nonbond_list(ia_2,is_2)%atom_type_number - !Convert epsilon to atomic units amu A^2/ps^2 - vdw_param1_table(itype_custom,jtype_custom) = kboltz * String_To_Double(line_array(3)) - vdw_param2_table(itype_custom,jtype_custom) = String_To_Double(line_array(4)) - !line_nbr = line_nbr + 1 - WRITE(logunit,'(A6,5x,A6,2x,T20,f10.4,T50,f10.4)') & - atom_type_list(itype_custom), atom_type_list(jtype_custom), & - vdw_param1_table(itype_custom,jtype_custom), vdw_param2_table(itype_custom,jtype_custom) - - END DO - END DO - END DO - END DO - RETURN - END IF - END DO - END IF - ENDDO + IF (int_vdw_style(1) == vdw_lj) THEN + ! Report parameters to logfile. Format is specific to vdw type. Add others here if + ! other than LJ potential is usd. + + WRITE(logunit,'(A6,5x,A6,2x,T20,f10.4,T50,f10.4)') & + atom_type_list(itype), atom_type_list(jtype), & + vdw_param1_table(itype,jtype), vdw_param2_table(itype,jtype) + + ENDIF + + IF (int_vdw_style(1) == born_cut_tail) THEN + ! Report parameters to logfile. Format is specific to vdw type. Add others here if + ! other than LJ potential is used. + + WRITE(logunit,'(A6,5x,A6,2x,T20,f18.6,T50,f10.4,T50,f18.6)') & + atom_type_list(itype), atom_type_list(jtype), & + vdw_param1_table(itype,jtype), vdw_param2_table(itype,jtype), vdw_param3_table(itype,jtype) !, vdw_param4_table(itype,jtype) !& + ! vdw_param5_table(itype,jtype), vdw_param6_table(itype,jtype) + + ENDIF + + + ENDDO - ENDDO - + ENDDO + WRITE(logunit,*) WRITE(logunit,*) '*** Completed construction of VDW interaction table ***' WRITE(logunit,*) diff --git a/Src/deletion.f90 b/Src/deletion.f90 index fc997146..2b02c444 100644 --- a/Src/deletion.f90 +++ b/Src/deletion.f90 @@ -267,7 +267,7 @@ SUBROUTINE Deletion(this_box,mcstep,randno) ! 4.4) Ewald energies - IF ( (int_charge_sum_style(this_box) == charge_ewald) .AND. & + IF ( (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian) .AND. & (has_charge(is)) ) THEN CALL Compute_Ewald_Reciprocal_Energy_Difference(alive,alive,is,this_box, & @@ -281,7 +281,7 @@ SUBROUTINE Deletion(this_box,mcstep,randno) ! 4.5) Long-range energy correction - IF (int_vdw_sum_style(this_box) == vdw_cut_tail) THEN + IF (int_vdw_sum_style(this_box) == vdw_cut_tail .or. int_vdw_sum_style(this_box) == born_cut_tail) THEN ! subtract off beads for this species nbeads_out(:) = nint_beads(:,this_box) @@ -367,13 +367,13 @@ SUBROUTINE Deletion(this_box,mcstep,randno) energy(this_box)%inter_vdw = energy(this_box)%inter_vdw - E_inter_vdw energy(this_box)%inter_q = energy(this_box)%inter_q - E_inter_qq - IF ( int_charge_sum_style(this_box) == charge_ewald .AND. & + IF ( (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian) .AND. & has_charge(is)) THEN energy(this_box)%ewald_reciprocal = E_reciprocal_move energy(this_box)%ewald_self = energy(this_box)%ewald_self + E_self_move END IF - IF ( int_vdw_sum_style(this_box) == vdw_cut_tail) THEN + IF ( int_vdw_sum_style(this_box) == vdw_cut_tail .or. int_vdw_sum_style(this_box) == born_cut_tail) THEN energy(this_box)%lrc = E_lrc END IF @@ -402,7 +402,7 @@ SUBROUTINE Deletion(this_box,mcstep,randno) ! CALL System_Energy_Check(1,mcstep,randno) ELSE - IF ( (int_charge_sum_style(this_box) == charge_ewald) .AND. & + IF ( (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian) .AND. & (has_charge(is)) ) THEN ! Restore cos_sum and sin_sum. Note that these were changed when ! difference in reciprocal energies was computed @@ -410,7 +410,7 @@ SUBROUTINE Deletion(this_box,mcstep,randno) sin_sum(:,this_box) = sin_sum_old(:,this_box) END IF - IF ( int_vdw_sum_style(this_box) == vdw_cut_tail ) THEN + IF ( int_vdw_sum_style(this_box) == vdw_cut_tail .or. int_vdw_sum_style(this_box) == born_cut_tail ) THEN ! Restore the total number of bead types nint_beads(:,this_box) = nbeads_out(:) END IF diff --git a/Src/dep.mk b/Src/dep.mk new file mode 100644 index 00000000..794a9408 --- /dev/null +++ b/Src/dep.mk @@ -0,0 +1,64 @@ +accept_or_reject.o: type_definitions.o run_variables.o random_generators.o +accumulate.o: run_variables.o simulation_properties.o +angle_distortion.o: type_definitions.o run_variables.o angle_dist_pick.o random_generators.o simulation_properties.o energy_routines.o pair_nrg_routines.o +angle_dist_pick.o: type_definitions.o run_variables.o file_names.o type_definitions.o run_variables.o random_generators.o file_names.o type_definitions.o run_variables.o file_names.o +atom_displacement.o: run_variables.o random_generators.o energy_routines.o simulation_properties.o +atoms_to_place.o: type_definitions.o run_variables.o file_names.o +chempot.o: run_variables.o energy_routines.o io_utilities.o random_generators.o rotation_routines.o type_definitions.o simulation_properties.o fragment_growth.o +clean_abort.o: file_names.o +compute_cell_dimensions.o: type_definitions.o run_variables.o io_utilities.o +create_intra_exclusion_table.o: run_variables.o type_definitions.o file_names.o +create_nonbond_table.o: run_variables.o type_definitions.o io_utilities.o file_names.o +cutNgrow.o: run_variables.o energy_routines.o fragment_growth.o random_generators.o simulation_properties.o io_utilities.o pair_nrg_routines.o +deletion.o: type_definitions.o run_variables.o random_generators.o simulation_properties.o energy_routines.o fragment_growth.o io_utilities.o +energy_routines.o: type_definitions.o run_variables.o file_names.o pair_nrg_routines.o random_generators.o run_variables.o run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o run_variables.o type_definitions.o run_variables.o type_definitions.o +file_names.o: +fragment_growth.o: run_variables.o energy_routines.o random_generators.o read_write_checkpoint.o rotation_routines.o io_utilities.o file_names.o energy_routines.o rotation_routines.o io_utilities.o file_names.o energy_routines.o run_variables.o random_generators.o run_variables.o random_generators.o io_utilities.o type_definitions.o rotation_routines.o +gcmc_control.o: io_utilities.o run_variables.o type_definitions.o file_names.o input_routines.o atoms_to_place.o angle_dist_pick.o energy_routines.o +gcmc_driver.o: run_variables.o random_generators.o file_names.o energy_routines.o read_write_checkpoint.o +gemc_control.o: io_utilities.o run_variables.o type_definitions.o file_names.o input_routines.o atoms_to_place.o angle_dist_pick.o energy_routines.o +gemc_driver.o: run_variables.o random_generators.o read_write_checkpoint.o energy_routines.o file_names.o io_utilities.o +gemc_nvt_volume.o: run_variables.o random_generators.o volume.o io_utilities.o energy_routines.o pair_nrg_routines.o +gemc_particle_transfer.o: run_variables.o random_generators.o simulation_properties.o energy_routines.o io_utilities.o fragment_growth.o pair_nrg_routines.o run_variables.o rotation_routines.o file_names.o +get_com.o: type_definitions.o run_variables.o type_definitions.o run_variables.o run_variables.o +get_internal_coords.o: type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o +grow_molecules.o: run_variables.o type_definitions.o random_generators.o file_names.o energy_routines.o fragment_growth.o simulation_properties.o io_utilities.o run_variables.o type_definitions.o random_generators.o file_names.o energy_routines.o fragment_growth.o +initialize.o: run_variables.o run_variables.o +input_routines.o: run_variables.o io_utilities.o file_names.o type_definitions.o random_generators.o energy_routines.o rotation_routines.o random_generators.o run_variables.o file_names.o file_names.o +insertion.o: run_variables.o energy_routines.o io_utilities.o random_generators.o rotation_routines.o fragment_growth.o +io_utilities.o: run_variables.o type_definitions.o +main.o: run_variables.o file_names.o io_utilities.o input_routines.o read_write_checkpoint.o energy_routines.o simulation_properties.o fragment_growth.o +mcf_control.o: io_utilities.o run_variables.o type_definitions.o input_routines.o atoms_to_place.o +min.o: run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o +minimum_image_separation.o: run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o +mpm_translate_rotate.o: type_definitions.o run_variables.o random_generators.o simulation_properties.o energy_routines.o pair_nrg_routines.o type_definitions.o run_variables.o random_generators.o simulation_properties.o energy_routines.o file_names.o pair_nrg_routines.o type_definitions.o random_generators.o +nptmc_control.o: io_utilities.o run_variables.o type_definitions.o file_names.o input_routines.o atoms_to_place.o angle_dist_pick.o energy_routines.o +nptmc_driver.o: run_variables.o random_generators.o file_names.o energy_routines.o read_write_checkpoint.o simulation_properties.o +nvtmc_control.o: io_utilities.o run_variables.o type_definitions.o file_names.o input_routines.o atoms_to_place.o angle_dist_pick.o energy_routines.o +nvtmc_driver.o: run_variables.o random_generators.o file_names.o energy_routines.o read_write_checkpoint.o +nvt_mc_fragment_control.o: io_utilities.o run_variables.o input_routines.o +nvt_mc_fragment_driver.o: run_variables.o random_generators.o energy_routines.o file_names.o run_variables.o random_generators.o +nvt_mc_ring_fragment.o: run_variables.o random_generators.o energy_routines.o run_variables.o random_generators.o fragment_growth.o io_utilities.o energy_routines.o run_variables.o random_generators.o energy_routines.o +pair_nrg_routines.o: type_definitions.o run_variables.o pair_nrg_variables.o +pair_nrg_variables.o: type_definitions.o +participation.o: type_definitions.o run_variables.o io_utilities.o file_names.o random_generators.o +precalculate.o: run_variables.o file_names.o energy_routines.o run_variables.o energy_routines.o type_definitions.o run_variables.o energy_routines.o +random_generators.o: type_definitions.o file_names.o run_variables.o type_definitions.o run_variables.o +read_write_checkpoint.o: run_variables.o file_names.o simulation_properties.o random_generators.o energy_routines.o +rigid_dihedral_change.o: type_definitions.o run_variables.o random_generators.o simulation_properties.o file_names.o energy_routines.o +rigid_insertion.o: type_definitions.o run_variables.o +rotate.o: type_definitions.o run_variables.o random_generators.o simulation_properties.o energy_routines.o file_names.o pair_nrg_routines.o +rotation_routines.o: run_variables.o random_generators.o +run_variables.o: type_definitions.o +save_revert_coordinates.o: type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o type_definitions.o run_variables.o +shell_gemc_nvt_volume.o: run_variables.o random_generators.o volume.o io_utilities.o energy_routines.o pair_nrg_routines.o +shell_gemc_particle_transfer.o: run_variables.o random_generators.o simulation_properties.o energy_routines.o io_utilities.o fragment_growth.o pair_nrg_routines.o volume.o +shell_relax.o: type_definitions.o run_variables.o random_generators.o simulation_properties.o energy_routines.o +simulation_properties.o: type_definitions.o run_variables.o +temp.o: +translate.o: type_definitions.o run_variables.o random_generators.o simulation_properties.o energy_routines.o pair_nrg_routines.o +type_definitions.o: +volume_change.o: type_definitions.o run_variables.o file_names.o random_generators.o energy_routines.o io_utilities.o +volume.o: run_variables.o type_definitions.o +write_properties.o: run_variables.o file_names.o energy_routines.o simulation_properties.o run_variables.o simulation_properties.o file_names.o +zig_by_omega.o: run_variables.o type_definitions.o random_generators.o file_names.o energy_routines.o fragment_growth.o diff --git a/Src/energy_routines.f90 b/Src/energy_routines.f90 index 0b31e70c..09da1ee7 100644 --- a/Src/energy_routines.f90 +++ b/Src/energy_routines.f90 @@ -1,3 +1,5 @@ + + !******************************************************************************** ! Cassandra - An open source atomistic Monte Carlo software package ! developed at the University of Notre Dame. @@ -968,6 +970,8 @@ SUBROUTINE Compute_Atom_Nonbond_Energy(this_atom,this_molecule,this_species,E_in vdw_in = vdw_cut_switch ELSEIF (int_vdw_sum_style(this_box) == vdw_mie) THEN vdw_in = vdw_mie + ELSEIF (int_vdw_sum_style(this_box) == born_cut_tail) THEN + vdw_in = born_cut_tail ENDIF IF (int_vdw_sum_style(this_box) /= vdw_in) THEN @@ -1093,12 +1097,24 @@ SUBROUTINE Compute_Atom_Nonbond_Energy(this_atom,this_molecule,this_species,E_in eps = eps * vdw_intra_scale(ia,this_atom,is) qsc = charge_intra_scale(ia,this_atom,is) ELSE + rij = SQRT(rijsq) + IF(vdw_in == born_cut_tail) THEN + SigoverR6 = rijsq*rijsq*rijsq + + IF(rij > vdw_param4_table(itype,jtype)) THEN + Eij_vdw = vdw_param1_table(itype,jtype)*DEXP(-vdw_param2_table(itype,jtype) * rij) - vdw_param3_table(itype,jtype)/sigoverR6 + ELSE + Eij_vdw = vdw_param1_table(itype,jtype) * 10000000.0_DP + END IF + + ELSE SigOverRsq = (sig**2) / rijsq SigOverR6 = SigOverRsq * SigOverRsq * SigOverRsq SigOverR12 = SigOverR6 * SigOverR6 Eij_vdw = 4.0_DP * eps * (SigOverR12 - SigOverR6) - + END IF + x = alpha_ewald(this_box) * rij T = 1.0_DP / (1.0_DP + P*x) xsq = x*x @@ -1517,6 +1533,10 @@ SUBROUTINE Pair_Energy & REAL(DP) :: roffsq_rijsq, roffsq_rijsq_sq, factor2, fscale REAL(DP) :: SigOverR, SigOverRn, SigOverRm, mie_coeff, rij, mie_n, mie_m, rij_shift, SigOverR_shift, SigOverRn_shift, SigOverRm_shift, rcut_vdw ! REAL(DP) :: Eij_vdw_check + + REAL(DP) :: A_eps,B_eps_r,C_eps_r6,r_born_max + REAL(DP) :: sqrt_r,r6 + Real(DP) :: qi,qj, qsc REAL(DP) :: this_lambda, RsqOverSig, R6OverSig, factorLJ REAL(DP) :: RsqOverSig_Shift, RsqOverSig6_Shift, factorLJ_Shift @@ -1629,6 +1649,33 @@ SUBROUTINE Pair_Energy & ENDIF LJ_12_6_calculation + + + Born_exp6_calculation: IF (int_vdw_style(1) == born_cut_tail) THEN + + sqrt_r = sqrt(rijsq) + R6 = rijsq*rijsq*rijsq + + !test buckingham exp-6 model alpha = 12 for water + A_eps = vdw_param1_table(itype,jtype) + B_eps_r = vdw_param2_table(itype,jtype) + C_eps_r6 = vdw_param3_table(itype,jtype) + r_born_max = vdw_param4_table(itype,jtype) + + IF (is == js .AND. im == jm) THEN + ! This controls 1-2, 1-3, and 1-4 interactions + A_eps = A_eps * vdw_intra_scale(ia,ja,is) + C_eps_r6 = C_eps_r6 * vdw_intra_scale(ia,ja,is) + ENDIF + + IF(rijsq > r_born_max*r_born_max) THEN + Eij_vdw = A_eps*EXP(-B_eps_r*sqrt_r) - C_eps_r6/R6 + ELSE + Eij_vdw = A_eps*100000000.0_DP + END IF + + ENDIF Born_exp6_calculation + ENDIF VDW_calculation @@ -1645,7 +1692,7 @@ SUBROUTINE Pair_Energy & qsc = charge_intra_scale(ia,ja,is) END IF Eij_qq = qsc*charge_factor*(qi*qj)/SQRT(rijsq) - ELSEIF (int_charge_sum_style(ibox) == charge_ewald .AND. ( .NOT. igas_flag) ) THEN + ELSEIF ((int_charge_sum_style(ibox) == charge_ewald .or. int_charge_sum_style(ibox) == charge_gaussian) .AND. ( .NOT. igas_flag) ) THEN ! Real space Ewald part this_box = molecule_list(im,is)%which_box CALL Ewald_Real(ia,im,is,qi,ja,jm,js,qj,rijsq,Eij_qq,this_box) @@ -1679,6 +1726,8 @@ SUBROUTINE Ewald_Real(ia,im,is,qi,ja,jm,js,qj,rijsq,Eij,ibox) INTEGER :: ia,im,is,ja,jm,js,ibox REAL(DP) :: qi,qj,qsc,rijsq,rij,erf_val REAL(DP) :: Eij + REAL(DP) :: alpha_ij, alpha_ia, alpha_ja + REAL(DP) :: fun_erfc qsc = 1.0_DP ibox = molecule_list(im,is)%which_box @@ -1697,9 +1746,41 @@ SUBROUTINE Ewald_Real(ia,im,is,qi,ja,jm,js,qj,rijsq,Eij,ibox) ! Come back to this later. rij = SQRT(rijsq) + + IF(int_charge_sum_style(ibox) == charge_ewald) THEN ! May need to protect against very small rijsq erf_val = 1.0_DP - erfc(alpha_ewald(ibox) * rij) Eij = (qi*qj/rij)*(qsc - erf_val)*charge_factor + + + ELSE IF(int_charge_sum_style(ibox) == charge_gaussian) THEN + + IF (abs(qi) > 1d-5 .and. abs(qj) > 1d-5) THEN + + alpha_ij = alp_ij(ia,is,ja,js) + + erf_val = 1.0_DP - derfc(alpha_ewald(ibox) * rij) + + ! erf_val = 1.0_DP - fun_erfc(alpha_ewald(ibox) * rij) + + ! print *, derfc(alpha_ewald(ibox) * rij), fun_erfc(alpha_ewald(ibox) * rij), rij + ! erf(alpha_ij*rij) + + IF(is .NE. js .or. im .NE. jm) THEN + Eij = qi*qj/rij*(1.0_DP - erfc(alpha_ij*rij) - erf_val)*charge_factor + ELSE + Eij = (qi*qj/rij)*(qsc - erf_val)*charge_factor + IF(rij<1d-6) Eij = -1.0_DP*qi*qj*alpha_ewald(ibox)*2.0/rootPI*charge_factor + END IF + + ELSE + + Eij = 0.0 + + END IF + + END IF + ! IF(en_flag) THEN ! WRITE(60,"(4I4,2F8.5,F24.12)") ia, ja, jm, js, qi,qj, Eij ! END IF @@ -2404,7 +2485,7 @@ SUBROUTINE Compute_Total_System_Energy(this_box,intra_flag,overlap) CALL Compute_Molecule_Improper_Energy(this_im,is,v_molecule_improper) intra_overlap = .FALSE. - IF (int_charge_sum_style(this_box) == charge_ewald) THEN + IF (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian) THEN CALL Compute_Molecule_Nonbond_Intra_Energy(this_im,is,vlj_molecule_intra,vqq_molecule_intra,intra_overlap, & v_molecule_selfrf) ELSE @@ -2423,7 +2504,7 @@ SUBROUTINE Compute_Total_System_Energy(this_box,intra_flag,overlap) v_improper = v_improper + v_molecule_improper v_intra_vdw = v_intra_vdw + vlj_molecule_intra v_intra_qq = v_intra_qq + vqq_molecule_intra - IF (int_charge_sum_style(this_box) == charge_ewald) THEN + IF (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian) THEN v_selfrf = v_selfrf + v_molecule_selfrf END IF @@ -2441,7 +2522,7 @@ SUBROUTINE Compute_Total_System_Energy(this_box,intra_flag,overlap) energy(this_box)%improper = energy(this_box)%improper + v_improper energy(this_box)%intra_vdw = energy(this_box)%intra_vdw + v_intra_vdw energy(this_box)%intra_q = energy(this_box)%intra_q + v_intra_qq - IF (int_charge_sum_style(this_box) == charge_ewald) THEN + IF (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian) THEN energy(this_box)%erf_self = energy(this_box)%erf_self + v_selfrf END IF END DO @@ -2610,7 +2691,7 @@ SUBROUTINE Compute_Total_System_Energy(this_box,intra_flag,overlap) ! Compute the reciprocal and self energy terms of the electrostatic energies if flag for Ewald is set. - IF (int_charge_sum_style(this_box) == charge_ewald) THEN + IF (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian) THEN ! Ewald reciprocal energy -- Note that the call changes global V_ewald_reciprocal(this_box) @@ -2630,7 +2711,7 @@ SUBROUTINE Compute_Total_System_Energy(this_box,intra_flag,overlap) END IF ! Long range correction if it is required - IF (int_vdw_sum_style(this_box) == vdw_cut_tail) THEN + IF (int_vdw_sum_style(this_box) == vdw_cut_tail .or. int_vdw_sum_style(this_box) == born_cut_tail) THEN CALL Compute_LR_Correction(this_box,e_lrc) ! add to the correction to the total energy of the system energy(this_box)%lrc = e_lrc @@ -2752,9 +2833,12 @@ SUBROUTINE Compute_LR_Correction(this_box, e_lrc) REAL(DP) :: epsij, sigij, sigij2, sigij6, sigij12 REAL(DP) :: e_lrc_ia_ja - - e_lrc = 0.0_DP + REAL(DP) :: A_eps, B_eps, C_eps, rcutoff + e_lrc = 0.0_DP + + IF(int_vdw_style(this_box) == vdw_lj) THEN + DO ia = 1, nbr_atomtypes e_lrc_ia_ja = 0.0_DP @@ -2779,6 +2863,27 @@ SUBROUTINE Compute_LR_Correction(this_box, e_lrc) e_lrc = e_lrc + REAL( nint_beads(ia,this_box), DP ) * e_lrc_ia_ja END DO + ELSE + + rcutoff = rcut3(this_box)**(1.0/3.0) + DO ia = 1, nbr_atomtypes + e_lrc_ia_ja = 0.0_DP + DO ja = 1, nbr_atomtypes + A_EPS = vdw_param1_table(ia,ja) + B_EPS = vdw_param2_table(ia,ja) + C_EPS = vdw_param3_table(ia,ja) + IF (B_EPS > tiny_number) THEN + e_lrc_ia_ja = e_lrc_ia_ja + & + (A_eps/(B_eps*B_eps*B_eps)*dexp(-1.0*B_eps*rcutoff)*(B_eps*B_eps*rcutoff*rcutoff + 2.0*B_eps*rcutoff + 2.0) - C_EPS/rcut3(this_box)/3.0 ) * REAL( nint_beads(ja,this_box), DP ) + END IF + END DO + + e_lrc = e_lrc + REAL( nint_beads(ia,this_box), DP ) * e_lrc_ia_ja + END DO + + ENDIF + + e_lrc = 2.0_DP * PI * e_lrc/box_list(this_box)%volume END SUBROUTINE Compute_LR_Correction @@ -2858,7 +2963,8 @@ SUBROUTINE Energy_Test(rijsq,get_vdw,get_qq,this_box) ELSE get_vdw = .FALSE. ENDIF - ELSEIF (int_vdw_sum_style(this_box) == vdw_cut .OR. int_vdw_sum_style(this_box) & + + ELSEIF (int_vdw_sum_style(this_box) == vdw_cut .OR. int_vdw_sum_style(this_box) & == vdw_cut_shift .OR. int_vdw_sum_style(this_box) == vdw_cut_tail) THEN IF (rijsq <= rcut_vdwsq(this_box)) THEN @@ -2888,10 +2994,24 @@ SUBROUTINE Energy_Test(rijsq,get_vdw,get_qq,this_box) ELSE get_vdw = .FALSE. END IF + + END IF + + ELSEIF(int_vdw_style(this_box) == born_cut_tail) THEN - ENDIF - - + IF (CBMC_flag) THEN + IF (rijsq <= rcut_cbmcsq) THEN + get_vdw = .TRUE. + ELSE + get_vdw = .FALSE. + ENDIF + ELSE + IF (rijsq <= rcut_vdwsq(this_box)) THEN + get_vdw = .TRUE. + ELSE + get_vdw = .FALSE. + ENDIF + END IF ELSE err_msg = "" @@ -2905,7 +3025,8 @@ SUBROUTINE Energy_Test(rijsq,get_vdw,get_qq,this_box) get_qq = .FALSE. ELSEIF (int_charge_style(this_box) == charge_coul) THEN - IF (int_charge_sum_style(this_box) == charge_cut .OR. int_charge_sum_style(this_box) == charge_ewald) THEN + IF (int_charge_sum_style(this_box) == charge_cut .OR. int_charge_sum_style(this_box) == charge_ewald .or. & + & int_charge_sum_style(this_box) == charge_gaussian ) THEN IF(CBMC_flag) THEN IF (rijsq <= rcut_cbmcsq) THEN get_qq = .TRUE. @@ -3046,13 +3167,13 @@ SUBROUTINE Compute_Forces(this_box) END DO imLOOP3 END DO - IF (int_charge_sum_style(this_box) == charge_ewald) THEN + IF (int_charge_sum_style(this_box) == charge_ewald .OR. int_charge_sum_style(this_box) == charge_gaussian ) THEN CALL Compute_System_Ewald_Reciprocal_Force(this_box) END IF - IF (int_vdw_sum_style(this_box) == vdw_cut_tail) THEN + IF (int_vdw_sum_style(this_box) == vdw_cut_tail .OR. int_vdw_sum_style(this_box) == born_cut_tail ) THEN CALL Compute_LR_Force(this_box,w_lrc) virial(this_box)%lrc = w_lrc @@ -3185,7 +3306,8 @@ SUBROUTINE Pair_Force & REAL(DP) :: roffsq_rijsq, roffsq_rijsq_sq, factor2, fscale REAL(DP) :: qi,qj, qsc, erf_val, this_lambda REAL(DP) :: rij, ewald_constant, exp_const, Wij_self - + + REAL(DP) :: A_eps, B_eps_r, C_eps_r6, sqrt_r, r_born_max, R6 !------------------------------------------------------------------------------------------ Wij_vdw = 0.0_DP @@ -3293,7 +3415,37 @@ SUBROUTINE Pair_Force & ! Add other potential types here ENDIF LJ_12_6_calculation - + + Born_cut_calculation: IF (int_vdw_style(1) == born_cut_tail) THEN + + sqrt_r = sqrt(rijsq) + R6 = rijsq*rijsq*rijsq + !R8 = R6*rijsq + !test buckingham exp-6 model alpha = 12 for water + A_eps = vdw_param1_table(itype,jtype) + B_eps_r = vdw_param2_table(itype,jtype) + C_eps_r6 = vdw_param3_table(itype,jtype) + r_born_max = vdw_param4_table(itype,jtype) + + IF (is == js .AND. im == jm) THEN + ! This controls 1-2, 1-3, and 1-4 interactions + A_eps = A_eps * vdw_intra_scale(ia,ja,is) + C_eps_r6 = C_eps_r6*vdw_intra_scale(ia,ja,is) + ENDIF + + IF(rijsq > r_born_max*r_born_max) THEN + Wij_vdw = A_eps*EXP(B_eps_r*(-sqrt_r))*sqrt_r*B_eps_r - C_eps_r6/R6*6.0 + ELSE + Wij_vdw = A_eps*100000000.0_DP + if(A_eps > tiny_number) then + write(*,*) 'error virial' + STOP + END IF + + END IF + ENDIF Born_cut_calculation + + ENDIF VDW_calculation qq_calculation: IF (get_qq) THEN @@ -3389,21 +3541,25 @@ SUBROUTINE Compute_LR_Force(this_box, w_lrc) ! ! !************************************************************************************************************** - INTEGER, INTENT(IN) :: this_box REAL(DP), INTENT(OUT) :: w_lrc - + INTEGER :: ia, ja REAL(DP) :: epsij, sigij REAL(DP) :: w_lrc_ia_ja + REAL(DP) :: A_epsij, B_epsij_r, C_epsij_R6 - w_lrc = 0.0_DP + REAL(DP) :: rcutof + w_lrc = 0.0_DP + + + IF(int_vdw_style(1) == vdw_lj) THEN DO ia = 1, nbr_atomtypes w_lrc_ia_ja = 0.0_DP - + DO ja = 1, nbr_atomtypes epsij = vdw_param1_table(ia,ja) @@ -3411,7 +3567,7 @@ SUBROUTINE Compute_LR_Force(this_box, w_lrc) w_lrc_ia_ja = w_lrc_ia_ja + nint_beads(ja,this_box) * epsij * ((2.0_DP / 3.0_DP * & sigij**12 / rcut9(this_box)) - (sigij**6 / rcut3(this_box))) - + END DO w_lrc = w_lrc + nint_beads(ia,this_box) * w_lrc_ia_ja @@ -3420,8 +3576,29 @@ SUBROUTINE Compute_LR_Force(this_box, w_lrc) w_lrc = 16.0_DP / 3.0_DP * PI * w_lrc / box_list(this_box)%volume + ELSE + + rcutof = rcut3(this_box)**(1.0/3.0) + DO ia = 1, nbr_atomtypes + w_lrc_ia_ja = 0.0_DP + DO ja = 1, nbr_atomtypes + A_EPSij = vdw_param1_table(ia,ja) + B_EPSij_R = vdw_param2_table(ia,ja) + C_EPSij_R6 = vdw_param3_table(ia,ja) + IF(A_epsij > tiny_number .and. B_EPSij_R>tiny_number) THEN + w_lrc_ia_ja = w_lrc_ia_ja + real(nint_beads(ja,this_box)) * & + (-C_epsij_r6/rcut3(this_box)/3.0) + END IF + END DO + w_lrc = w_lrc + real(nint_beads(ia,this_box)) * w_lrc_ia_ja + END DO + w_lrc = 2.0_DP / 3.0_DP * PI * w_lrc / box_list(this_box)%volume + END IF + END SUBROUTINE Compute_LR_Force + + SUBROUTINE Compute_System_Ewald_Reciprocal_Force(this_box) !***************************************************************************************** ! This subroutine computes the long range forces due to electrostatics @@ -4027,5 +4204,591 @@ SUBROUTINE Compute_System_Ewald_Reciprocal_Energy2(this_box) END SUBROUTINE Compute_System_Ewald_Reciprocal_Energy2 + + + + ! polarizable force field related calculation + + !==================Force and Torque Calculation=================================================== + + + SUBROUTINE compute_molecule_inter_overlap(im, is, this_box, overlap) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: im, is, this_box + LOGICAL, INTENT(OUT) :: overlap + + INTEGER :: ispecies, imolecule, this_locate, ja, ia + REAL(DP) :: rxijp, ryijp, rzijp, rxij, ryij, rzij, rijsq + LOGICAL :: shared_overlap + + overlap = .FALSE. + shared_overlap =.FALSE. + speciesLoop: DO ispecies = 1, nspecies + + + moleculeLoop: DO imolecule = 1, nmolecules(ispecies) + + this_locate = locate(imolecule,ispecies) + + IF (ispecies == is .AND. this_locate == im) CYCLE moleculeLOOP + IF (molecule_list(this_locate,ispecies)%which_box /= this_box) CYCLE moleculeLOOP + IF (.NOT. molecule_list(this_locate,ispecies)%live) CYCLE moleculeLOOP + IF (shared_overlap) CYCLE moleculeLOOP + + DO ia = 1, natoms(is) + DO ja = 1, natoms(ispecies) + + rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,this_locate,ispecies)%rxp + ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,this_locate,ispecies)%ryp + rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,this_locate,ispecies)%rzp + CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp,rxij,ryij,rzij) + + rijsq = rxij*rxij + ryij*ryij + rzij*rzij + + IF(rijsq < rcut_lowsq) THEN + overlap = .TRUE. + END IF + + END DO + END DO + + END DO moleculeLoop + + IF (overlap) shared_overlap = .TRUE. + END DO speciesLoop + + END SUBROUTINE compute_molecule_inter_overlap + + + + + SUBROUTINE Compute_Molecule_Nonbond_Inter_Force(im, is, F1, F2, F3, E_inter_vdw, E_inter_qq, overlap, enecal) + + IMPLICIT NONE + + INTEGER, INTENT(IN):: im, is + REAL(DP), INTENT(OUT) :: F1(maxval(natoms)), F2(maxval(natoms)), F3(maxval(natoms)) + REAL(DP), INTENT(OUT) :: E_inter_vdw, E_inter_qq + LOGICAL, INTENT(IN) :: enecal + LOGICAL :: overlap + + + !--------------------------------------------------------------------------------------------------- + + INTEGER :: this_box, ispecies, imolecule, this_locate, j + + REAL(DP) :: Eij_vdw, Eij_qq + REAL(DP) :: eps + REAL(DP) :: rcom, rx, ry, rz + + LOGICAL :: get_interaction + + INTEGER :: locate_1, locate_2 + + LOGICAL :: l_pair_store + LOGICAL :: my_overlap, shared_overlap + REAL(DP) :: F_x(maxval(natoms)), F_y(maxval(natoms)), F_z(maxval(natoms)) + + + E_inter_vdw = 0.0 + E_inter_qq = 0.0 + DO j = 1, maxval(natoms) + F1(j) = 0.0_DP; F_x(j) = 0.0 + F2(j) = 0.0_DP; F_y(j) = 0.0 + F3(j) = 0.0_DP; F_z(j) = 0.0 + END DO + + + overlap = .FALSE. + my_overlap = .FALSE. + shared_overlap = .false. + + this_box = molecule_list(im,is)%which_box + + speciesLoop: DO ispecies = 1, nspecies + + + moleculeLoop: DO imolecule = 1, nmolecules(ispecies) + + IF(shared_overlap) CYCLE + + this_locate = locate(imolecule,ispecies) + + IF (ispecies == is .AND. this_locate == im) CYCLE moleculeLOOP + + IF (molecule_list(this_locate,ispecies)%which_box /= this_box) CYCLE moleculeLOOP + + IF(.NOT. molecule_list(this_locate,ispecies)%live) CYCLE moleculeLOOP + + CALL Check_Interaction(im,is,this_locate,ispecies,get_interaction,rcom,rx,ry,rz) + + IF (.NOT. get_interaction) CYCLE moleculeLOOP + + CALL Compute_Molecule_Pair_Interaction_Force(im,is,this_locate,ispecies,this_box, & + F_x, F_y, F_z, Eij_vdw, Eij_qq, my_overlap, enecal) + + DO j = 1, maxval(natoms) + F1(j) = F1(j) + F_x(j) + F2(j) = F2(j) + F_y(j) + F3(j) = F3(j) + F_z(j) + END DO + + IF(my_overlap) shared_overlap = .TRUE. + + E_inter_vdw = E_inter_vdw + Eij_vdw + E_inter_qq = E_inter_qq + Eij_qq + + + END DO moleculeLoop + + IF(shared_overlap) THEN + overlap = .TRUE. + RETURN + END IF + + END DO speciesLoop + + END SUBROUTINE Compute_Molecule_Nonbond_Inter_Force + + SUBROUTINE Compute_Molecule_Pair_Interaction_Force(im_1,is_1,im_2,is_2,this_box, F_x, F_y, F_z, vlj_pair, vqq_pair, overlap, enecal) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: im_1, is_1, im_2, is_2, this_box + LOGICAL, INTENT(IN) :: enecal + REAL(DP), INTENT(OUT) :: F_x(maxval(natoms)), F_y(maxval(natoms)), F_z(maxval(natoms)), vlj_pair, vqq_pair + LOGICAL, INTENT(OUT) :: overlap + !---------------------------------------------------------------------------------------------------------- + + INTEGER :: ia, ja + + REAL(DP) :: rxijp, ryijp, rzijp, rxij, ryij, rzij, rijsq, rij, Eij_vdw, Eij_qq + + LOGICAL :: get_vdw, get_qq + + INTEGER :: locate_im_1, locate_im_2, old_style + + REAL(DP) :: eps, sig, sigoverrsq,sigOverR6, SigOverR12 + REAL(DP) :: qi, qj, qsc, r_rc + REAL(DP) :: p_rc, p_rc_prime, du_dr, R_C + REAL(DP) :: A_eps, B_eps_r, C_eps_r6, D_eps_r8, sqrt_r, r6, r8, r_born_max + REAL(DP) :: alpha_ij, erf_val, fun_erfc + INTEGER :: itype, jtype + INTEGER :: this_box1 + + F_x = 0.0_DP + F_y = 0.0_DP + F_z = 0.0_DP + vlj_pair = 0.0_DP + vqq_pair = 0.0_DP + qsc = 1.0_DP + + DO ia = 1, natoms(is_1) + + IF (.NOT. atom_list(ia,im_1,is_1)%exist) CYCLE + + DO ja = 1, natoms(is_2) + ! Obtain the minimum image separation + + IF ( .NOT. atom_list(ja,im_2,is_2)%exist) CYCLE + + rxijp = atom_list(ia,im_1,is_1)%rxp - atom_list(ja,im_2,is_2)%rxp + ryijp = atom_list(ia,im_1,is_1)%ryp - atom_list(ja,im_2,is_2)%ryp + rzijp = atom_list(ia,im_1,is_1)%rzp - atom_list(ja,im_2,is_2)%rzp + + ! 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 + + IF(rijsq < rcut_lowsq) THEN + overlap = .TRUE. + RETURN + END IF + + CALL Energy_Test(rijsq,get_vdw,get_qq,this_box) + + + IF (get_vdw) THEN + + IF (int_vdw_style(1) == born_cut_tail) THEN + + IF (atom_list(ia,im_1,is_1)%exist .AND. atom_list(ja,im_2,is_2)%exist) THEN + + itype = nonbond_list(ia,is_1)%atom_type_number + jtype = nonbond_list(ja,is_2)%atom_type_number + rij = sqrt(rijsq) + R6 = rijsq*rijsq*rijsq + + A_eps = vdw_param1_table(itype,jtype) + B_eps_r = vdw_param2_table(itype,jtype) + C_eps_r6 = vdw_param3_table(itype,jtype) + r_born_max = vdw_param4_table(itype,jtype) + + IF (is_1 == is_2 .AND. im_1 == im_2) THEN + ! This controls 1-2, 1-3, and 1-4 interactions + A_eps = A_eps * vdw_intra_scale(ia,ja,is_1) + C_eps_r6 = C_eps_r6* vdw_intra_scale(ia,ja,is_1) + ENDIF + + IF(rij > r_born_max) THEN + du_dr = (B_eps_r*A_eps*dexp(-B_eps_r*rij)*rij - C_eps_r6/R6)/rijsq + F_x(ia) = F_x(ia) + du_dr*rxij + F_y(ia) = F_y(ia) + du_dr*ryij + F_z(ia) = F_z(ia) + du_dr*rzij + ELSE + F_x(ia) = F_x(ia) + A_eps*100000000.0_DP + F_y(ia) = F_y(ia) + A_eps*100000000.0_DP + F_z(ia) = F_z(ia) + A_eps*100000000.0_DP + END IF + + + IF (enecal) THEN + Eij_vdw = A_eps*dexp(-B_eps_r*rij) - C_eps_r6/R6 + + if(rij < r_born_max) then + Eij_vdw = A_eps*1d7; + end if + + vlj_pair = Vlj_pair + Eij_vdw + + END IF + + END IF + + END IF + + IF (int_vdw_style(1) == vdw_lj) THEN + + IF (atom_list(ia,im_1,is_1)%exist .AND. atom_list(ja,im_2,is_2)%exist) THEN + + itype = nonbond_list(ia,is_1)%atom_type_number + jtype = nonbond_list(ja,is_2)%atom_type_number + + eps = vdw_param1_table(itype,jtype) + sig = vdw_param2_table(itype,jtype) + + IF (is_1 == is_2 .AND. im_1 == im_2) THEN + eps = eps * vdw_intra_scale(ia,ja,is_1) + ENDIF + + SigOverRsq = sig*sig/rijsq + SigOverR6 = SigOverRsq * SigOverRsq * SigOverRsq + SigOverR12 = SigOverR6 * SigOverR6 + + du_dr = 4.0_DP*eps*(2.0*SigOverR12 - SigOverR6)/Rijsq + + F_x(ia) = F_x(ia) + du_dr*rxij + F_y(ia) = F_y(ia) + du_dr*ryij + F_z(ia) = F_z(ia) + du_dr*rzij + IF (enecal) THEN + Eij_vdw = 4.0_DP*eps*(SigOverR12 - SigOverR6) + vlj_pair = Vlj_pair + Eij_vdw + END IF + END IF + END IF + END IF + + IF (get_qq) THEN + qsc = 1.0_DP + IF (is_1 == is_2 .AND. im_1 == im_2) THEN + ! This controls 1-2, 1-3, and 1-4 interactions + qsc = charge_intra_scale(ia,ja,is_1) + ENDIF + qi = nonbond_list(ia,is_1)%charge + qj = nonbond_list(ja,is_2)%charge + IF (abs(qi) > 1D-5 .and. abs(qj) > 1D-5) THEN + + rij = sqrt(rijsq) + du_dr = qsc/rij*qi*qj + + F_x(ia) = F_x(ia) + rxij*reaf(int(rij/reaf_width),this_box)*du_dr ! use the generated reaction field table for quick force calculation + F_y(ia) = F_y(ia) + ryij*reaf(int(rij/reaf_width),this_box)*du_dr + F_z(ia) = F_z(ia) + rzij*reaf(int(rij/reaf_width),this_box)*du_dr + + IF(int_charge_sum_style(this_box) == charge_gaussian) THEN + + alpha_ij = alp_ij(ia,is_1,ja,is_2) + + erf_val = 1.0_DP - erfc(alpha_ij*rij) + DU_Dr = -1.0*qsc*charge_factor*qi*qj*( -(erf_val-1.0)/rijsq + 2.0/rootPI*alpha_ij*exp(-alpha_ij*alpha_ij*rijsq)/rij) + F_x(ia) = F_x(ia) + DU_DR*rxij/rij + F_y(ia) = F_y(ia) + DU_DR*ryij/rij + F_z(ia) = F_z(ia) + DU_DR*rzij/rij + END IF + + END IF + this_box1 = molecule_list(im_2,is_2)%which_box + IF(this_box1 .NE. this_box) THEN + print *, 'incorrect box index in pair force calculation' + STOP + END IF + IF (enecal) THEN + CALL Ewald_Real(ia,im_1,is_1,qi,ja,im_2,is_2,qj,rijsq,Eij_qq,this_box) + vqq_pair = vqq_pair + Eij_qq + END IF + + END IF + + END DO + + END DO + + END SUBROUTINE Compute_Molecule_Pair_Interaction_Force + + + + SUBROUTINE Compute_Molecule_Torque(im, is, this_box, F_x, F_y, F_z, M_xx, M_yy, M_zz) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: im, is, this_box + REAL(DP), INTENT(IN) :: F_x(maxval(natoms)), F_y(maxval(natoms)), F_z(maxval(natoms)) + REAL(DP), INTENT(OUT) :: M_xx, M_yy, M_zz + REAL(DP) :: rxijp, ryijp, rzijp, rxij, ryij, rzij + INTEGER :: i, ia + + M_xx = 0.0_DP + M_yy = 0.0_DP + M_zz = 0.0_DP + + DO ia = 1, natoms(is) + + rxijp = atom_list(ia,im,is)%rxp - molecule_list(im,is)%xcom + ryijp = atom_list(ia,im,is)%ryp - molecule_list(im,is)%ycom + rzijp = atom_list(ia,im,is)%rzp - molecule_list(im,is)%zcom + + CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp,rxij,ryij,rzij) + + M_zz = -F_x(ia)*ryij + F_y(ia)*rxij + M_zz + M_yy = -F_z(ia)*rxij + F_x(ia)*rzij + M_yy + M_xx = -F_y(ia)*rzij + F_z(ia)*ryij + M_xx + + END DO + + END SUBROUTINE Compute_Molecule_Torque + + + + SUBROUTINE Compute_drude_cos_sin_sum(this_box) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: this_box + + INTEGER :: i, is, im, ia, this_locate + REAL(DP) :: hdotr + REAL(DP) :: k_sq_i, ek_beta + REAL(DP) :: cos_val, sin_val + + !$OMP PARALLEL WORKSHARE DEFAULT(SHARED) + cos_sum(:,this_box) = 0.0_DP + sin_sum(:,this_box) = 0.0_DP + !$OMP END PARALLEL WORKSHARE + + DO is = 1, nspecies + + IF ( .NOT. has_charge(is)) CYCLE + + DO im = 1, nmolecules(is) + + this_locate = locate(im,is) + + IF( .NOT. molecule_list(this_locate,is)%live) CYCLE + IF( molecule_list(this_locate,is)%which_box /= this_box ) CYCLE + + !$OMP PARALLEL DO DEFAULT(SHARED)& + !$OMP SCHEDULE(STATIC)& + !$OMP PRIVATE(i,ia,hdotr, k_sq_i, cos_val, sin_val, ek_beta) + + DO i = 1, nvecs(this_box) + + cos_mol(i,this_locate) = 0.0_DP + sin_mol(i,this_locate) = 0.0_DP + + DO ia = 1, natoms(is) + + IF(abs(nonbond_list(ia,is)%charge)<1D-5) CYCLE !this site does not have charge + ! compute hdotr + hdotr = hx(i,this_box) * atom_list(ia,this_locate,is)%rxp + & + hy(i,this_box) * atom_list(ia,this_locate,is)%ryp + & + hz(i,this_box) * atom_list(ia,this_locate,is)%rzp + + cos_mol(i,this_locate) = cos_mol(i,this_locate) + nonbond_list(ia,is)%charge * DCOS(hdotr) + sin_mol(i,this_locate) = sin_mol(i,this_locate) + nonbond_list(ia,is)%charge * DSIN(hdotr) + + END DO + + cos_sum(i,this_box) = cos_sum(i,this_box) + cos_mol(i,this_locate) + sin_sum(i,this_box) = sin_sum(i,this_box) + sin_mol(i,this_locate) + END DO + + !$OMP END PARALLEL DO + END DO + + END DO + + END SUBROUTINE Compute_drude_cos_sin_sum + + + SUBROUTINE Compute_molecule_drude_force(is1,im1,ia1,this_box, Fxd, Fyd, Fzd) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: is1, im1, ia1, this_box + REAL(DP), INTENT(OUT):: Fxd, Fyd, Fzd + + INTEGER :: i, is, im, ia, this_locate, k + + REAL(DP) :: un, const_val, this_scaling + REAL(DP) :: charge_i, charge_j + REAL(DP) :: ewald_constant, exp_const + REAL(DP) :: rxijp, ryijp, rzijp, rxij, ryij, rzij, rijsq, rij + REAL(DP) :: erf_val, forced, hdotr, qsc + REAL(DP) :: alpha_ij, exp_factor + REAL(DP) :: force_factor, exp1 + REAL(DP) :: fun_erfc + + Fxd = 0.0_DP + Fyd = 0.0_DP + Fzd = 0.0_DP + charge_i = nonbond_list(ia1,is1)%charge + + DO i = 1, nvecs(this_box) + hdotr = hx(i,this_box) * atom_list(ia1,im1,is1)%rxp + & + hy(i,this_box) * atom_list(ia1,im1,is1)%ryp + & + hz(i,this_box) * atom_list(ia1,im1,is1)%rzp + + hdotr = DSIN(hdotr)*cos_sum(i,this_box) - DCOS(hdotr)*sin_sum(i,this_box) + hdotr = hdotr*2.0*Cn(i,this_box) + + Fxd = Fxd + hdotr*hx(i,this_box) + Fyd = Fyd + hdotr*hy(i,this_box) + Fzd = Fzd + hdotr*hz(i,this_box) + END DO + + DO is = 1, nspecies + + IF ( .NOT. has_charge(is)) CYCLE + + DO im = 1, nmolecules(is) + + this_locate = locate(im,is) + + IF( .NOT. molecule_list(this_locate,is)%live) CYCLE + IF( molecule_list(this_locate,is)%which_box /= this_box ) CYCLE + + DO ia = 1, natoms(is) + + IF(abs(nonbond_list(ia1,is1)%charge) < 1d-5 .or. abs(nonbond_list(ia,is)%charge) < 1d-5 ) CYCLE + + rxijp = atom_list(ia1,im1,is1)%rxp - atom_list(ia,this_locate,is)%rxp + ryijp = atom_list(ia1,im1,is1)%ryp - atom_list(ia,this_locate,is)%ryp + rzijp = atom_list(ia1,im1,is1)%rzp - atom_list(ia,this_locate,is)%rzp + + CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp,rxij,ryij,rzij) + rijsq = rxij*rxij + ryij*ryij + rzij*rzij + + + IF (rijsq < rcut_coulsq(this_box) ) THEN + + rij = DSQRT(rijsq) + + IF(int_charge_sum_style(this_box) == charge_ewald) THEN + ewald_constant = 2.0_DP * alpha_ewald(this_box) / rootPI + exp_const = DEXP(-1.0_DP*alpha_ewald(this_box)*alpha_ewald(this_box)*rijsq) + ! May need to protect against very small rijsq + erf_val = 1.0_DP - erfc(alpha_ewald(this_box) * rij) + qsc = 1.0 + IF(is == is1 .and. this_locate==im1) then + qsc = charge_intra_scale(ia,ia1,is1) + END IF + forced = (ewald_constant*exp_const + (qsc - erf_val )/rij ) * nonbond_list(ia,is)%charge + IF (is .NE. is1 .or. this_locate .NE. im1 .or. ia .NE. ia1) THEN + IF(rij > 1d-6) THEN + Fxd = Fxd + forced*rxij/rijsq + Fyd = Fyd + forced*ryij/rijsq + Fzd = Fzd + forced*rzij/rijsq + END IF + END IF + + ELSE IF(int_charge_sum_style(this_box) == charge_gaussian) THEN + + ewald_constant = 2.0_DP * alpha_ewald(this_box) / rootPI + exp_const = DEXP(-1.0_DP*alpha_ewald(this_box)*alpha_ewald(this_box)*rijsq) + + erf_val = 1.0_DP - derfc(alpha_ewald(this_box) * rij) + + IF (is == is1 .and. this_locate == im1) THEN + if(rij>1d-5 .and. ia .NE. ia1) then + qsc = charge_intra_scale(ia,ia1,is1) + force_factor = (ewald_constant*exp_const + (qsc - erf_val )/rij ) * nonbond_list(ia,is)%charge / rijsq + Fxd = Fxd + force_factor*rxij + Fyd = Fyd + force_factor*ryij + Fzd = Fzd + force_factor*rzij + end if + ELSE + alpha_ij = alp_ij(ia,is,ia1,is1) + exp1 = exp(-1.0_DP*alpha_ij*alpha_ij*rijsq) + force_factor = (1.0_DP - derfc(alpha_ij*rij))/rij - erf_val/rij + ewald_constant*exp_const - 2.0*alpha_ij*exp1/rootPI + force_factor = force_factor/rijsq*nonbond_list(ia,is)%charge + Fxd = Fxd + force_factor*rxij + Fyd = Fyd + force_factor*ryij + Fzd = Fzd + force_factor*rzij + END IF + END IF + END IF + END DO + END DO + END DO + + Fxd = Fxd * charge_i * charge_factor + Fyd = Fyd * charge_i * charge_factor + Fzd = Fzd * charge_i * charge_factor + + + if(ia1 - 1 .LE. 0) then + print *, 'drude particle should be attached to the previous atom, fatal error' + end if + + rxijp = atom_list(ia1-1,im1,is1)%rxp - atom_list(ia1,im1,is1)%rxp + ryijp = atom_list(ia1-1,im1,is1)%ryp - atom_list(ia1,im1,is1)%ryp + rzijp = atom_list(ia1-1,im1,is1)%rzp - atom_list(ia1,im1,is1)%rzp + CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp,rxij,ryij,rzij) + Fxd = Fxd + rxij*nonbond_list(ia1,is1)%pol_alpha + Fyd = Fyd + ryij*nonbond_list(ia1,is1)%pol_alpha + Fzd = Fzd + rzij*nonbond_list(ia1,is1)%pol_alpha + + END SUBROUTINE Compute_molecule_drude_force + END MODULE Energy_Routines +REAL FUNCTION fun_erfc(x) + + USE Type_Definitions + + IMPLICIT NONE + + REAL(DP) :: x + + Real(DP), Parameter :: a1 = 0.254829592 + Real(DP), Parameter :: a2 = -0.284496736 + Real(DP), Parameter :: a3 = 1.421413741 + Real(DP), Parameter :: a4 = -1.453152027 + Real(DP), Parameter :: a5 = 1.061405429 + Real(DP), Parameter :: pp = 0.3275911 + + Real(DP) :: tt, exp_val + + tt = 1.0d0/(1.0d0+pp*x) + exp_val = exp(-x*x) + fun_erfc = tt*(a1+tt*(a2+tt*(a3+tt*(a4+tt*a5))))*exp_val + + RETURN + +END FUNCTION fun_erfc diff --git a/Src/gemc_driver.f90 b/Src/gemc_driver.f90 index 3288b3df..b5f0b258 100644 --- a/Src/gemc_driver.f90 +++ b/Src/gemc_driver.f90 @@ -87,10 +87,9 @@ SUBROUTINE GEMC_Driver i = i + 1 ! We will select a move from Golden Sampling scheme - rand_no = rranf() which_step = i - + IF (rand_no <= cut_trans) THEN IF(.NOT. openmp_flag) THEN @@ -99,7 +98,11 @@ SUBROUTINE GEMC_Driver !$ time_s = omp_get_wtime() END IF - CALL Translate(this_box,which_step) + IF(shell_mpm) THEN + CALL Translate_MP(this_box,which_step) + ELSE + CALL Translate(this_box,which_step) + END IF IF(.NOT. openmp_flag) THEN CALL cpu_time(time_e) @@ -108,7 +111,7 @@ SUBROUTINE GEMC_Driver END IF movetime(imove_trans) = movetime(imove_trans) + time_e - time_s - + ELSE IF ( rand_no <= cut_rot) THEN IF(.NOT. openmp_flag) THEN @@ -116,8 +119,12 @@ SUBROUTINE GEMC_Driver ELSE !$ time_s = omp_get_wtime() END IF - - CALL Rotate(this_box) + + IF(shell_mpm) THEN + CALL Rotate_MP(this_box) + ELSE + CALL Rotate(this_box) + END IF IF(.NOT. openmp_flag) THEN CALL cpu_time(time_e) @@ -159,8 +166,12 @@ SUBROUTINE GEMC_Driver ELSE + IF(shell_mpm) THEN + CALL Shell_GEMC_NVT_Volume(this_box, other_box) + ELSE CALL GEMC_NVT_Volume(this_box, other_box) volume_move = .TRUE. + END IF END IF @@ -198,7 +209,11 @@ SUBROUTINE GEMC_Driver !$ time_s = omp_get_wtime() END IF - CALL GEMC_Particle_Transfer(this_box,other_box) + IF(shell_mpm) THEN + CALL Shell_GEMC_Particle_Transfer(this_box,other_box) + ELSE + CALL GEMC_Particle_Transfer(this_box,other_box) + END IF IF(.NOT. openmp_flag) THEN CALL cpu_time(time_e) @@ -248,6 +263,7 @@ SUBROUTINE GEMC_Driver END IF + ! Accumulate averages. Note that the averages must ! be updated for both the boxes for the moves that ! affect both the boxes. These moves are volume @@ -302,15 +318,14 @@ SUBROUTINE GEMC_Driver END IF END IF - IF(write_flag) THEN + IF(write_flag) THEN DO ibox = 1, nbr_boxes CALL Write_Properties(i,ibox) CALL Reset(ibox) - DO j = 1,nspecies IF(cpcollect) WRITE(*,*) 'Chem Pot spec:',j,'=',chpot(j,ibox) / ntrials(j,ibox)%cpcalc @@ -319,7 +334,8 @@ SUBROUTINE GEMC_Driver END DO END IF - + + write_flag = .FALSE. IF(.NOT. timed_run) THEN diff --git a/Src/input_routines.f90 b/Src/input_routines.f90 index 57b74ec0..9c7dc489 100644 --- a/Src/input_routines.f90 +++ b/Src/input_routines.f90 @@ -427,7 +427,7 @@ SUBROUTINE Get_Pair_Style vdw_style(ibox) = line_array(1) WRITE(logunit,'(A,2x,A,A,I3)') ' VDW style used is: ',vdw_style(ibox), 'in box:', ibox - IF (vdw_style(ibox) /= 'NONE') THEN + IF (vdw_style(ibox) /= 'NONE' .AND. vdw_style(ibox) /= 'Born' ) THEN int_vdw_style(ibox) = vdw_lj vdw_sum_style(ibox) = line_array(2) WRITE(logunit,'(A,2x,A,A,I3)') ' VDW sum style is: ',vdw_sum_style(ibox), 'in box:', ibox @@ -509,14 +509,16 @@ SUBROUTINE Get_Pair_Style ELSEIF (vdw_sum_style(ibox) == 'mie') THEN int_vdw_sum_style(ibox) = vdw_mie - rcut_vdw(ibox) = String_To_Double(line_array(3)) + rcut_vdw(ibox) = String_To_Double(line_array(3)) WRITE(logunit,'(A,2x,F7.3, A)') ' rcut = ',rcut_vdw(ibox), ' Angstrom' WRITE(logunit,'(A)') 'Mie potential used for VDW' - + ELSE err_msg(1) = 'Improper specification of vdw_sum_style' CALL Clean_Abort(err_msg,'Get_Pairstyle') ENDIF + + IF (rcut_vdw(ibox) > MIN(box_list(ibox)%face_distance(1)/2.0_DP, & box_list(ibox)%face_distance(2)/2.0_DP, box_list(ibox)%face_distance(3)/2.0_DP)) THEN @@ -528,7 +530,15 @@ SUBROUTINE Get_Pair_Style ENDIF - + ELSEIF (vdw_style(ibox) == 'Born') THEN !The Born potential (exp-6 potential) + + int_vdw_style(ibox) = born_cut_tail + vdw_sum_style(ibox) = line_array(2) + int_vdw_sum_style(ibox) = born_cut_tail + rcut_vdw(ibox) = String_To_Double(line_array(3)) + rcut3(ibox) = rcut_vdw(ibox) * rcut_vdw(ibox) * rcut_vdw(ibox) + rcut9(ibox) = rcut3(ibox) * rcut3(ibox) * rcut3(ibox) + !The Born potential is cut and tail correction will be added ELSE @@ -579,8 +589,14 @@ SUBROUTINE Get_Pair_Style WRITE(logunit,'(A,2x,F7.3, A)') ' rcut = ',rcut_coul(ibox), ' Angstrom' - ELSEIF (charge_sum_style(ibox) == 'Ewald') THEN - int_charge_sum_style(ibox) = charge_ewald + ELSEIF (charge_sum_style(ibox) == 'Ewald' .or. charge_sum_style(ibox) == 'Gaussian_charge' ) THEN + + IF(charge_sum_style(ibox) == 'Ewald') THEN + int_charge_sum_style(ibox) = charge_ewald + ELSE + int_charge_sum_style(ibox) = charge_gaussian + END IF + rcut_coul(ibox) = String_To_Double(line_array(3)) IF (l_half_len_cutoff(ibox)) THEN @@ -747,6 +763,8 @@ SUBROUTINE Get_Mixing_Rules WRITE(logunit,'(A)') 'Lorentz-Berthelot mixing rule specified' ELSEIF (mix_rule == 'geometric') THEN WRITE(logunit,'(A)') 'Geometric mixing rule specified' + ELSE IF (mix_rule == 'kong') THEN + WRITE(logunit,'(A)') 'kong mixing rule' ELSEIF (mix_rule == 'custom') THEN WRITE(logunit,'(A)') 'Custom mixing rule specified' ELSE @@ -1115,7 +1133,8 @@ SUBROUTINE Get_Molecule_Info IF ( .NOT. (int_sim_type == sim_frag .OR. int_sim_type == sim_ring )) THEN CALL Get_Insertion_Style(is) END IF - CALL Get_Atom_Info(is) + + CALL Get_Atom_Info(is) CALL Get_Bond_Info(is) CALL Get_Angle_Info(is) CALL Get_Dihedral_Info(is) @@ -1482,6 +1501,24 @@ SUBROUTINE Get_Atom_Info(is) ! Set number of vdw parameters nbr_vdw_params = 2 + + ELSE IF (nonbond_list(ia,is)%vdw_potential_type == 'Born') THEN !Born (exp-6) potential + + nonbond_list(ia,is)%vdw_param(1) = String_To_Double(line_array(7)) ! A + nonbond_list(ia,is)%vdw_param(2) = String_To_Double(line_array(8)) ! B + nonbond_list(ia,is)%vdw_param(3) = String_To_Double(line_array(9)) ! C + nonbond_list(ia,is)%vdw_param(4) = String_To_Double(line_array(10)) ! rmax + WRITE(logunit,'(A,T25,F18.4)') ' Epsilon / kB in K:', & + nonbond_list(ia,is)%vdw_param(1) + WRITE(logunit,'(A,T25,F10.4)') ' Sigma B / A^-1:', & + nonbond_list(ia,is)%vdw_param(2) + WRITE(logunit,'(A,T25,F18.4)') ' EpsilonC / kB in K:', & + nonbond_list(ia,is)%vdw_param(3) + + nonbond_list(ia,is)%vdw_param(1) = kboltz* nonbond_list(ia,is)%vdw_param(1) + nonbond_list(ia,is)%vdw_param(3) = kboltz* nonbond_list(ia,is)%vdw_param(3) + nbr_vdw_params = 4 + ELSEIF (nonbond_list(ia,is)%vdw_potential_type == 'NONE') THEN WRITE(logunit,'(A,I6,1x,I6)') & 'No VDW potential assigned to atom, species: ',ia,is @@ -1510,6 +1547,24 @@ SUBROUTINE Get_Atom_Info(is) nexo_atoms(is) = nexo_atoms(is) + 1 exo_atom_ids(nexo_atoms(is),is) = ia END IF + + + ! Input the width of the Gaussian Charge, input the polarizability of atom (polarizability of core) (the last 2 entry) + IF ( (nbr_entries .GE. 9 .and. (.NOT. nonbond_list(ia,is)%ring_atom) .and. nonbond_list(ia,is)%vdw_potential_type .NE. 'Born') & + & .or. (nbr_entries .GT. 10 .and. nonbond_list(ia,is)%vdw_potential_type == 'Born') )THEN + + nonbond_list(ia,is)%g_alpha = String_To_Double(line_array(nbr_entries-1)) + + nonbond_list(ia,is)%g_alpha = 1.0/nonbond_list(ia,is)%g_alpha/sqrt(2.0) + + nonbond_list(ia,is)%pol_alpha = String_To_Double(line_array(nbr_entries)) + + IF(nonbond_list(ia,is)%pol_alpha > tiny_number) nonbond_list(ia,is)%pol_alpha = 1389.32198*100.0*nonbond_list(ia,is)%charge*nonbond_list(ia,is)%charge/nonbond_list(ia,is)%pol_alpha + + WRITE(logunit,*) 'gaussian width of atom ', ia, 'species', is, ' equal to ', 1.0/sqrt(2.0)/nonbond_list(ia,is)%g_alpha, String_To_Double(line_array(nbr_entries-1)) + WRITE(logunit,*) 'polarizability of atom ', ia, 'species', is, ' equal to ', String_To_Double(line_array(nbr_entries)) + + END IF WRITE(logunit,*) @@ -1530,6 +1585,7 @@ SUBROUTINE Get_Atom_Info(is) ENDIF ENDDO + WRITE(logunit,'(A, T40, I4,A, T45, I4)') 'Total number of ring atoms in species', is, ' is', nring_atoms(is) WRITE(logunit,'(A, T40, I4,A, T45, I4)') 'Total number of exo atoms in species', is, ' is', nexo_atoms(is) @@ -1635,6 +1691,19 @@ SUBROUTINE Get_Bond_Info(is) ! Set number of bond parameters nbr_bond_params = 1 + ELSE IF (bond_list(ib,is)%bond_potential_type == 'harmonic' ) THEN + bond_list(ib,is)%int_bond_type = int_harmonic + WRITE(logunit,'(A,I6,1x,I6, A, I4)') & + 'Harmonic bond between atoms: ',bond_list(ib,is)%atom1, bond_list(ib,is)%atom2, & + 'in species', is + ! Harmonic bond + bond_list(ib,is)%bond_param(1) = String_To_Double(line_array(5)) + bond_list(ib,is)%bond_param(2) = String_To_Double(line_array(6)) + WRITE(logunit,*) 'Fixed bond length, in A:',bond_list(ib,is)%bond_param(1) + + ! Set number of bond parameters + nbr_bond_params = 2 + ELSE err_msg = "" err_msg(1) = 'bond_potential type improperly specified in mcf file' @@ -3451,6 +3520,10 @@ SUBROUTINE Get_Box_Info ALLOCATE(P_inst(nbr_boxes),P_ideal(nbr_boxes)) + ALLOCATE(P_inst_plus(nbr_boxes),P_inst_minus(nbr_boxes)) + + constant_vol = .FALSE. + W_tensor_charge(:,:,:) = 0.0_DP W_tensor_recip(:,:,:) = 0.0_DP W_tensor_vdw(:,:,:) = 0.0_DP diff --git a/Src/insertion.f90 b/Src/insertion.f90 index 429e8bc4..9075f989 100644 --- a/Src/insertion.f90 +++ b/Src/insertion.f90 @@ -379,7 +379,7 @@ SUBROUTINE Insertion(this_box,mcstep,randno) ! 3.5) Ewald energies - IF ( (int_charge_sum_style(this_box) == charge_ewald) .AND. & + IF ( (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian ) .AND. & (has_charge(is)) ) THEN CALL Compute_Ewald_Reciprocal_Energy_Difference(alive,alive,is,this_box, & @@ -394,7 +394,7 @@ SUBROUTINE Insertion(this_box,mcstep,randno) ! 3.6) Long-range energy correction - IF (int_vdw_sum_style(this_box) == vdw_cut_tail) THEN + IF (int_vdw_sum_style(this_box) == vdw_cut_tail .or. int_vdw_sum_style(this_box) == born_cut_tail ) THEN ! increase number of integer beads nbeads_in = nint_beads(:,this_box) @@ -476,13 +476,13 @@ SUBROUTINE Insertion(this_box,mcstep,randno) energy(this_box)%inter_vdw = energy(this_box)%inter_vdw + E_inter_vdw energy(this_box)%inter_q = energy(this_box)%inter_q + E_inter_qq - IF ( int_charge_sum_style(this_box) == charge_ewald .AND. & + IF ( (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian) .AND. & has_charge(is)) THEN energy(this_box)%ewald_reciprocal = E_reciprocal_move energy(this_box)%ewald_self = energy(this_box)%ewald_self + E_self_move END IF - IF (int_vdw_sum_style(this_box) == vdw_cut_tail) THEN + IF (int_vdw_sum_style(this_box) == vdw_cut_tail .or. int_vdw_sum_style(this_box) == born_cut_tail) THEN energy(this_box)%lrc = e_lrc END IF @@ -495,7 +495,7 @@ SUBROUTINE Insertion(this_box,mcstep,randno) atom_list(:,alive,is)%exist = .FALSE. molecule_list(alive,is)%molecule_type = int_none - IF ( int_charge_sum_style(this_box) == charge_ewald .AND. & + IF ( (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian) .AND. & has_charge(is) ) THEN ! Restore cos_sum and sin_sum. Note that these were changed when the ! difference in reciprocal energies was computed. @@ -503,7 +503,7 @@ SUBROUTINE Insertion(this_box,mcstep,randno) sin_sum(:,this_box) = sin_sum_old(:,this_box) END IF - IF ( int_vdw_sum_style(this_box) == vdw_cut_tail ) THEN + IF ( int_vdw_sum_style(this_box) == vdw_cut_tail .or. int_vdw_sum_style(this_box) == born_cut_tail ) THEN ! Restore the total number of bead types nint_beads(:,this_box) = nbeads_in(:) END IF diff --git a/Src/main.f90 b/Src/main.f90 index 7142f941..f2d865eb 100644 --- a/Src/main.f90 +++ b/Src/main.f90 @@ -89,7 +89,7 @@ PROGRAM Main ! !$ include 'omp_lib.h' INTEGER(4) :: count - INTEGER :: i, j, is, im, ia, this_im, ibox, nmol_is, this_box, int_phi + INTEGER :: i, j, is, is1, im, ia, ia1, this_im, ibox, nmol_is, this_box, int_phi INTEGER :: alive, t_num INTEGER :: nyears, nmonths, ndays, nhours, nmin, nsec, nms @@ -97,7 +97,7 @@ PROGRAM Main CHARACTER(120) filename1 CHARACTER(80) :: name - LOGICAL :: overlap, cbmc_overlap, superbad, inside_ch + LOGICAL :: overlap, cbmc_overlap, superbad, inside_ch, conv_drude REAL(DP) :: attempt_prob, phi REAL(DP) :: E_st_vdw, E_st_qq, W_st_vdw, W_st_qq, e_lrc, w_lrc @@ -261,6 +261,30 @@ PROGRAM Main CALL Initialize(i) END DO + IF(int_charge_sum_style(1) == charge_gaussian) THEN + ALLOCATE(alp_ij(MAXVAL(natoms), nspecies, MAXVAL(natoms), nspecies), alp_ia(MAXVAL(natoms), nspecies, nbr_boxes)) + alp_ij = 0.0_DP + alp_ia = 0.0_DP + DO is = 1, nspecies + DO ia = 1, natoms(is) + + DO is1 = 1, nspecies + DO ia1 = 1, natoms(is1) + alp_ij(ia,is,ia1,is1) = nonbond_list(ia,is)%g_alpha * nonbond_list(ia1,is1)%g_alpha & + & / dsqrt(nonbond_list(ia,is)%g_alpha * nonbond_list(ia,is)%g_alpha + nonbond_list(ia1,is1)%g_alpha * nonbond_list(ia1,is1)%g_alpha) + END DO + END DO + + DO ibox = 1, nbr_boxes + alp_ia(ia,is,ibox) = nonbond_list(ia,is)%g_alpha * alpha_ewald(ibox) & + & / dsqrt(nonbond_list(ia,is)%g_alpha * nonbond_list(ia,is)%g_alpha + alpha_ewald(ibox) * alpha_ewald(ibox)) + END DO + + END DO + END DO + END IF + + WRITE(*,*) 'Beginning Cassandra Simulation' WRITE(*,*) @@ -290,7 +314,7 @@ PROGRAM Main ! Ewald stuff - IF ( int_charge_sum_style(1) == charge_ewald) THEN + IF ( int_charge_sum_style(1) == charge_ewald .or. int_charge_sum_style(1) == charge_gaussian ) THEN ALLOCATE(hx(maxk,nbr_boxes),hy(maxk,nbr_boxes),hz(maxk,nbr_boxes), & hsq(maxk,nbr_boxes), Cn(maxk,nbr_boxes),Stat=AllocateStatus) @@ -413,14 +437,18 @@ PROGRAM Main ! compute total system energy overlap = .FALSE. - + CALL Initialize_MP + DO i = 1, nbr_boxes IF ( start_type == 'make_config' .and. & - int_vdw_sum_style(i) == vdw_cut_tail) & + int_vdw_sum_style(i) == vdw_cut_tail .or. int_vdw_sum_style(this_box) == born_cut_tail) & CALL Compute_Beads(i) + IF(shell_mpm) CALL SHELL_RELAX(i, conv_drude) CALL Compute_Total_System_Energy(i,.TRUE.,overlap) - + print *, 'multiparticle move with shell particle', shell_mpm + print *, 'EM for box',i,'convergence = ', conv_drude + IF (overlap) THEN ! overlap was detected between two atoms so abort the program @@ -429,7 +457,12 @@ PROGRAM Main err_msg(2) = 'Start type '//start_type CALL Clean_Abort(err_msg,'Main') END IF + END DO + + + + ! write out the initial energy components to the log file @@ -449,19 +482,19 @@ PROGRAM Main WRITE(logunit,'(X,A,T30,F20.3)') 'Intra nonbond vdw is', energy(i)%intra_vdw WRITE(logunit,'(X,A,T30,F20.3)') 'Intra nonbond elec is', energy(i)%intra_q WRITE(logunit,'(X,A,T30,F20.3)') 'Inter molecule vdw is', energy(i)%inter_vdw - IF (int_vdw_sum_style(i) == vdw_cut_tail) THEN + IF (int_vdw_sum_style(i) == vdw_cut_tail .or. int_vdw_sum_style(i) == born_cut_tail) THEN WRITE(logunit,'(X,A,T30,F20.3)') 'Long range correction is', energy(i)%lrc END IF WRITE(logunit,'(X,A,T30,F20.3)') 'Inter molecule q is', energy(i)%inter_q WRITE(logunit,'(X,A,T30,F20.3)') 'Reciprocal ewald is', energy(i)%ewald_reciprocal WRITE(logunit,'(X,A,T30,F20.3)') 'Self ewald is', energy(i)%ewald_self - IF( int_charge_sum_style(i) == charge_ewald) WRITE(logunit,'(X,A,T30,I20)') 'Number of vectors is', nvecs(i) + IF( int_charge_sum_style(i) == charge_ewald .or. int_charge_sum_style(i) == charge_gaussian) WRITE(logunit,'(X,A,T30,I20)') 'Number of vectors is', nvecs(i) WRITE(logunit,'(X,A59)') '***********************************************************' WRITE(logunit,*) END DO - + ! obtain the number of molecules of each species in a simulation box ! only if it isn't a new configuration. diff --git a/Src/mpm_translate_rotate.f90 b/Src/mpm_translate_rotate.f90 new file mode 100644 index 00000000..9d4a7614 --- /dev/null +++ b/Src/mpm_translate_rotate.f90 @@ -0,0 +1,730 @@ +!******************************************************************************** +! Cassandra - An open source atomistic Monte Carlo software package +! developed at the University of Notre Dame. +! http://cassandra.nd.edu +! Prof. Edward Maginn +! Copyright (2013) University of Notre Dame du Lac +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see . +!******************************************************************************** + + + + +!******************************************************************************** +! This subroutine performs translate/rotate move for all the particles in this_box. +! This is a force biased move. +SUBROUTINE Translate_MP(this_box,mc_step) +!******************************************************************************** + + USE Type_Definitions + USE Run_Variables + USE Random_Generators + USE Simulation_Properties + USE Energy_Routines + USE Pair_Nrg_Routines + + IMPLICIT NONE + INTEGER, INTENT(IN) :: this_box + INTEGER :: is, im, i, j, nmolecules_species, alive, total_mols, ibox, nmols_box, this_species + INTEGER :: mc_step, ia + + REAL(DP) :: delta_e, p_acc + REAL(DP) :: E_intra, E_vdw, E_qq + REAL(DP) :: E_vdw_move, E_qq_move + REAL(DP) :: E_intra_qq, E_intra_vdw + REAL(DP) :: E_intra_qq_move, E_intra_vdw_move + REAL(DP) :: E_intra_qq_alive, E_intra_vdw_alive, E_intra_qq_move_alive, E_intra_vdw_move_alive + REAL(DP) :: E_bond, E_bond_move, E_bond_alive + REAL(DP) :: E_reciprocal_move + REAL(DP) :: Energy_old_reciprocal + REAL(DP) :: E_vdw_move_alive,E_qq_move_alive, E_vdw_alive, E_qq_alive + + REAL(DP) :: success_ratio + REAL(DP) :: rcut_small + + LOGICAL :: inter_overlap, overlap, total_overlap, intra_overlap + LOGICAL :: energy_cal + LOGICAL :: conv_drude + + REAL(DP) :: attempt_p + + REAL(DP) :: gamma_x_a , gamma_y_a , gamma_z_a + REAL(DP) :: S_new_old, S_old_new, S_is + REAL(kind=8) :: g_random + REAL(kind=8) :: dev_b, mean_a, zero, one + REAL(kind=8) :: dispp_max + +! REAL(dp) :: g_random +! REAL(dp) :: dev_b, mean_a +! REAL(dp) :: dispp_max + + E_vdw_move = 0.0_DP + E_qq_move = 0.0_DP + E_vdw = 0.0_DP + E_qq = 0.0_DP + E_reciprocal_move = 0.0_DP + attempt_p = 1.0_DP + + mean_a = 0.D0 + one = 1.D0 + + total_overlap = .FALSE. + + conv_Drude = .TRUE. + + + + total_mols = SUM(nmols(:,this_box)) + + ! If there are no molecules in the box then return to driver + IF (total_mols == 0) RETURN + + + rcut_small = MIN(rcut_vdw(this_box),rcut_coul(this_box)) + IF (rcut_small < 1.0) rcut_small = rcut_vdw(this_box) + + ! update the trial counter for this molecule + + tot_trials(this_box) = tot_trials(this_box) + 1 + DO is = 1, nspecies + ntrials(is,this_box)%displacement = ntrials(is,this_box)%displacement + 1 + END DO + + E_vdw = 0.0; E_qq = 0.0; inter_overlap = .FALSE. + + energy_cal =.FALSE. + + DO is = 1, nspecies + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP SCHEDULE(STATIC) & + !$OMP PRIVATE(im, alive) + DO im = 1, nmolecules(is) !nmols(is,this_box) + alive = locate(im,is) + if(molecule_list(alive,is)%which_box /= this_box) cycle + if(.NOT. molecule_list(alive,is)%live) cycle + CALL Save_Old_Cartesian_Coordinates(alive,is) + call Compute_Molecule_Nonbond_Inter_Force(alive, is, Fxx(is, alive, :), Fyy(is, alive, :), Fzz(is, alive, :), & + & E_vdw_alive, E_qq_alive, inter_overlap, energy_cal) + + IF (inter_overlap) THEN + WRITE(*,*) 'Error in the old configuration -- > overlap', im ,is + STOP + END IF + END DO + !$OMP END PARALLEL DO + END DO + + + E_vdw = energy(this_box)%inter_vdw + E_qq = energy(this_box)%inter_q + E_intra_qq = energy(this_box)%intra_q + E_intra_vdw = energy(this_box)%intra_vdw + E_bond = energy(this_box)%bond + E_intra = energy(this_box)%intra + + ! Compute the displacement vector for all partciles + ! dispp_max is an adjustable parameter for multiparticle move + + S_new_old = 0.0; S_old_new = 0.0 + + LOOP1:DO is = 1, nspecies + S_is = 0.0 + dispp_max = max_disp(is,this_box) + dev_b = dsqrt(dispp_max*2.0_DP) + if(dispp_max*1000.0 < 1d-7) cycle + + LOOP2:DO im = 1, nmolecules(is) !nmols(is,this_box) + alive = locate(im,is) + if(molecule_list(alive,is)%which_box /= this_box) cycle + if(.NOT. molecule_list(alive,is)%live) cycle + gamma_x_a = g_random(mean_a, dev_b) + gamma_y_a = g_random(mean_a, dev_b) + gamma_z_a = g_random(mean_a, dev_b) + dmpmx(is,alive) = beta(this_box)*dispp_max*sum(Fxx(is, alive, :)) + gamma_x_a + dmpmy(is,alive) = beta(this_box)*dispp_max*sum(Fyy(is, alive, :)) + gamma_y_a + dmpmz(is,alive) = beta(this_box)*dispp_max*sum(Fzz(is, alive, :)) + gamma_z_a + atom_list(:,alive,is)%rxp = atom_list(:,alive,is)%rxp + dmpmx(is,alive) + atom_list(:,alive,is)%ryp = atom_list(:,alive,is)%ryp + dmpmy(is,alive) + atom_list(:,alive,is)%rzp = atom_list(:,alive,is)%rzp + dmpmz(is,alive) + molecule_list(alive,is)%xcom = molecule_list(alive,is)%xcom + dmpmx(is,alive) + molecule_list(alive,is)%ycom = molecule_list(alive,is)%ycom + dmpmy(is,alive) + molecule_list(alive,is)%zcom = molecule_list(alive,is)%zcom + dmpmz(is,alive) + S_is = S_is + (beta(this_box)*dispp_max*sum(Fxx(is, alive, :)) - dmpmx(is,alive))*& + (beta(this_box)*dispp_max*sum(Fxx(is, alive, :)) - dmpmx(is,alive)) + S_is = S_is + (beta(this_box)*dispp_max*sum(Fyy(is, alive, :)) - dmpmy(is,alive))*& + (beta(this_box)*dispp_max*sum(Fyy(is, alive, :)) - dmpmy(is,alive)) + S_is = S_is + (beta(this_box)*dispp_max*sum(Fzz(is, alive, :)) - dmpmz(is,alive))*& + (beta(this_box)*dispp_max*sum(Fzz(is, alive, : )) - dmpmz(is,alive)) + + END DO LOOP2 + S_old_new = S_old_new + S_is/dispp_max + END DO LOOP1 + S_old_new = S_old_new /4.D0*(-1.D0) + + + E_vdw_move = 0.0; E_qq_move = 0.0; inter_overlap = .FALSE.; total_overlap = .FALSE. + + LOOP4: DO is = 1, nspecies + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP SCHEDULE(STATIC) & + !$OMP PRIVATE(im, alive, inter_overlap) & + !$OMP REDUCTION(.OR.:total_overlap) + LOOP3: DO im = 1, nmolecules(is) !nmols(is,this_box) + alive = locate(im,is) + if(molecule_list(alive,is)%which_box /= this_box) cycle + if(.NOT. molecule_list(alive,is)%live) cycle + CALL Fold_Molecule(alive,is,this_box) + CALL compute_molecule_inter_overlap(alive,is,this_box,inter_overlap) + IF (inter_overlap) THEN + total_overlap = .TRUE. + END IF + END DO LOOP3 + !$OMP END PARALLEL DO + END DO LOOP4 + + + IF (total_overlap ) THEN + DO is = 1, nspecies + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP SCHEDULE(STATIC) & + !$OMP PRIVATE(im, alive) + DO im = 1, nmolecules(is) + alive = locate(im,is) + if(molecule_list(alive,is)%which_box /= this_box) cycle + if(.NOT. molecule_list(alive,is)%live) cycle + CALL Revert_Old_Cartesian_Coordinates(alive,is) + END DO + !$OMP END PARALLEL DO + END DO + END IF + + + IF (.NOT. total_overlap) THEN + + conv_drude = .TRUE. + IF(shell_mpm) CALL shell_relax(this_box, conv_drude) ! perform electrostatic energy minimization + + E_vdw_move = 0.0; E_qq_move = 0.0 + E_intra_vdw_move = 0.0; E_intra_qq_move = 0.0 + E_bond_move = 0.0 + + energy_cal = .TRUE. + IF (conv_drude) THEN + DO is = 1, nspecies + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP SCHEDULE(STATIC) & + !$OMP PRIVATE(im, alive) & + !$OMP PRIVATE(E_vdw_move_alive, E_qq_move_alive,inter_overlap) & + !$OMP PRIVATE(E_intra_vdw_move_alive, E_intra_qq_move_alive,E_bond_alive,intra_overlap) & + !$OMP REDUCTION(+:E_vdw_move,E_qq_move, E_intra_qq_move, E_intra_vdw_move, E_bond_move) + + DO im = 1, nmolecules(is) !nmols(is,this_box) + alive = locate(im,is) + if(molecule_list(alive,is)%which_box /= this_box) cycle + if(.NOT. molecule_list(alive,is)%live) cycle + CALL Compute_Molecule_Nonbond_Inter_Force(alive, is, Fxx(is, alive, :), Fyy(is, alive, :), Fzz(is, alive, :), & + & E_vdw_move_alive, E_qq_move_alive, inter_overlap, energy_cal) + CALL Compute_Molecule_Nonbond_Intra_Energy(alive,is,E_intra_vdw_move_alive,E_intra_qq_move_alive,intra_overlap) + CALL Compute_Molecule_Bond_Energy(alive,is,E_bond_alive) + E_vdw_move = E_vdw_move_alive + E_vdw_move + E_qq_move = E_qq_move_alive + E_qq_move + E_intra_qq_move = E_intra_qq_move_alive + E_intra_qq_move + E_intra_vdw_move = E_intra_vdw_move_alive + E_intra_vdw_move + E_bond_move = E_bond_move + E_bond_alive + END DO + + !$OMP END PARALLEL DO + + END DO + + LOOP22: DO is = 1, nspecies + S_is = 0.0 + dispp_max = max_disp(is,this_box) + if(dispp_max*1000.0 < 1d-7) cycle + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP SCHEDULE(STATIC) & + !$OMP PRIVATE(im, alive) & + !$OMP REDUCTION(+:S_is) + LOOP23: DO im = 1, nmolecules(is) !nmols(is,this_box) + alive = locate(im,is) + if(molecule_list(alive,is)%which_box /= this_box) cycle + if(.NOT. molecule_list(alive,is)%live) cycle + S_is = S_is + (beta(this_box)*dispp_max*sum(Fxx(is, alive, :)) + dmpmx(is,alive))*& + (beta(this_box)*dispp_max*sum(Fxx(is, alive, :)) + dmpmx(is,alive)) + S_is = S_is + (beta(this_box)*dispp_max*sum(Fyy(is, alive, :)) + dmpmy(is,alive))*& + (beta(this_box)*dispp_max*sum(Fyy(is, alive, :)) + dmpmy(is,alive)) + S_is = S_is + (beta(this_box)*dispp_max*sum(Fzz(is, alive, :)) + dmpmz(is,alive))*& + (beta(this_box)*dispp_max*sum(Fzz(is, alive, :)) + dmpmz(is,alive)) + END DO LOOP23 + !$OMP END PARALLEL DO + S_new_old = S_is/dispp_max + S_new_old + END DO LOOP22 + S_new_old = S_new_old / 4.0 * (-1.0) + + END IF + + END IF + + IF ((int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian) .and. .NOT. total_overlap .and. conv_drude) THEN + Energy_old_reciprocal = energy(this_box)%ewald_reciprocal + CALL Compute_System_Ewald_Reciprocal_Energy2(this_box) + END IF + + delta_e = ( E_vdw_move/2.0 - E_vdw ) + ( E_qq_move/2.0 - E_qq ) + delta_e = delta_e + E_bond_move - E_bond + E_intra_vdw_move + E_intra_qq_move - E_intra_qq - E_intra_vdw + delta_e = delta_e + energy(this_box)%ewald_reciprocal - Energy_old_reciprocal + + p_acc = MIN( 1.0_DP, ( DEXP(-beta(this_box) * delta_e + S_new_old - S_old_new ) ) ) + + IF ( rranf() <= p_acc .and. .NOT. total_overlap .and. conv_drude ) THEN + + energy(this_box)%inter_vdw = E_vdw_move/2.0 + energy(this_box)%inter_q = E_qq_move/2.0 + energy(this_box)%intra_q = E_intra_qq_move + energy(this_box)%intra_vdw = E_intra_vdw_move + energy(this_box)%bond = E_bond_move + energy(this_box)%intra = E_bond_move + energy(this_box)%total = energy(this_box)%total + delta_e + nsuccess(:,this_box)%displacement = nsuccess(:,this_box)%displacement + 1 + nequil_success(:,this_box)%displacement = nequil_success(:,this_box)%displacement + 1 + + ELSE + + IF (.NOT. total_overlap) THEN + DO is = 1, nspecies + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP SCHEDULE(STATIC) & + !$OMP PRIVATE(im, alive) + DO im = 1, nmolecules(is) !nmols(is,this_box) + alive = locate(im,is) + if(molecule_list(alive,is)%which_box /= this_box) cycle + if(.NOT. molecule_list(alive,is)%live) cycle + + CALL Revert_Old_Cartesian_Coordinates(alive,is) + END DO + !$OMP END PARALLEL DO + END DO + + IF (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian .and. conv_drude) THEN + Energy(this_box)%ewald_reciprocal = Energy_old_reciprocal + END IF + + END IF + + END IF + + + + DO is = 1, nspecies + IF ( MOD(ntrials(is,this_box)%displacement,nupdate) == 0 ) THEN + IF ( int_run_style == run_equil ) THEN + success_ratio = REAL(nequil_success(is,this_box)%displacement,DP)/REAL(nupdate,DP) + ELSE + success_ratio = REAL(nsuccess(is,this_box)%displacement,DP)/REAL(ntrials(is,this_box)%displacement,DP) + END IF + WRITE(logunit,*) + WRITE(logunit,'(A,2X,I3,2X,A,I2,2x,A,2X,F8.5)')' Multiparticle Translation ratio for species ', is , 'for box', this_box, 'is: ', success_ratio + IF ( int_run_style == run_equil ) THEN + ! check if the acceptace is close to 0.33 + nequil_success(is,this_box)%displacement = 0 + IF ( success_ratio < 0.001 ) THEN + max_disp(is,this_box) = 0.1_DP*max_disp(is,this_box) + ELSE + max_disp(is,this_box) = MIN(rcut_small, 3.0_DP*success_ratio*max_disp(is,this_box)) + END IF + IF(max_disp(is,this_box) .gt. 25.0) max_disp(is,this_box) = 25.0 + WRITE(logunit,'(A,1X,I1,1X,A,1X,I1)') 'Parameter A for multiparticle translation of species', is,' updated in box', this_box + WRITE(logunit,'(A,2X,F10.5)') 'is', max_disp(is,this_box) + WRITE(logunit,*) + END IF + END IF +END DO + + +END SUBROUTINE Translate_MP + + + + +SUBROUTINE Rotate_MP(this_box, mc_step) + + + USE Type_Definitions + USE Run_Variables + USE Random_Generators + USE Simulation_Properties + USE Energy_Routines + USE File_Names + USE Pair_Nrg_Routines + IMPLICIT NONE + +! !$ include 'omp_lib.h' + INTEGER, INTENT(IN) :: this_box + INTEGER :: mc_step + INTEGER :: is, im, ia, nmolecules_species, alive, i, nmols_box, total_mols, this_species + INTEGER :: ibox + + REAL(DP) :: E_intra, E_vdw, E_qq + REAL(DP) :: E_intra_qq, E_intra_vdw + REAL(DP) :: E_vdw_alive, E_qq_alive, E_vdw_move_alive, E_qq_move_alive, Energy_old_reciprocal + REAL(DP) :: E_intra_qq_alive, E_intra_vdw_alive, E_intra_vdw_move_alive, E_intra_qq_move_alive + REAL(DP) :: E_intra_qq_move, E_intra_vdw_move, E_vdw_move, E_qq_move + REAL(DP) :: E_bond, E_bond_move, E_bond_alive, E_reciprocal_move + REAL(DP) :: delta_e, p_acc, success_ratio + REAL(DP) :: attempt_p + LOGICAL :: inter_overlap + LOGICAL :: overlap,total_overlap, energy_cal,intra_overlap + LOGICAL :: conv_drude + + REAL(DP) :: S_new_old , S_old_new, S_is + REAL(DP) :: gamma_x_a, gamma_y_a, gamma_z_a + REAL(kind=8) :: g_random + REAL(kind=8) :: mean_a, dev_b + REAL(kind=8):: rott_max + + + attempt_p = 1.0_DP + + inter_overlap = .FALSE. + E_vdw_move = 0.0_DP + E_vdw = 0.0_DP + E_qq_move = 0.0_DP + E_qq = 0.0_DP + E_reciprocal_move = 0.0_DP + mean_a = 0.D0 + + total_overlap = .FALSE. + conv_Drude = .TRUE. + + total_mols = SUM(nmols(:,this_box)) + ! If there are no molecules in the box then return + + IF (total_mols == 0) RETURN + + tot_trials(this_box) = tot_trials(this_box) + 1 + DO is = 1, nspecies + ntrials(is,this_box)%rotation = ntrials(is,this_box)%rotation + 1 + END DO + + ! Evaluate force and torque on each molecule + + E_vdw = 0.0; E_qq = 0.0; inter_overlap = .FALSE. + + energy_cal = .FALSE. + + LOOP31: DO is = 1, nspecies + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP SCHEDULE(STATIC) & + !$OMP PRIVATE(im, alive) + LOOP32: DO im = 1, nmolecules(is) ! nmols(is,this_box) + alive = locate(im,is) + if(molecule_list(alive,is)%which_box /= this_box) cycle + if(.NOT. molecule_list(alive,is)%live) cycle + + CALL Save_Old_Cartesian_Coordinates(alive,is) + CALL Compute_Molecule_Nonbond_Inter_Force(alive, is, Fxx(is, alive, :), Fyy(is, alive, :), Fzz(is, alive, :), & + & E_vdw_alive, E_qq_alive, inter_overlap, energy_cal) + IF (inter_overlap) THEN + WRITE(*,*) 'Error in the old configuration "in rotation" -- > overlap' + STOP + END IF + CALL Compute_Molecule_Torque(alive, is, this_box, Fxx(is, alive,:), Fyy(is, alive,:), Fzz(is, alive,:), Mxx(is, alive), Myy(is, alive), Mzz(is,alive) ) + END DO LOOP32 + !$OMP END PARALLEL DO + END DO LOOP31 + + + E_vdw = energy(this_box)%inter_vdw + E_qq = energy(this_box)%inter_q + E_intra_qq = energy(this_box)%intra_q + E_intra_vdw = energy(this_box)%intra_vdw + E_bond = energy(this_box)%bond + E_intra = energy(this_box)%intra + + S_new_old = 0.0; S_old_new = 0.0 + + LOOP33: DO is = 1, nspecies + S_is = 0.0_DP + rott_max = max_rot(is,this_box) + dev_b = dsqrt(rott_max*2.0) +! if(is .gt. 2) cycle + LOOP34: DO im = 1, nmolecules(is) !nmols(is,this_box) + alive = locate(im,is) + if(molecule_list(alive,is)%which_box /= this_box) cycle + if(.NOT. molecule_list(alive,is)%live) cycle + gamma_x_a = g_random(mean_a, dev_b) + gamma_y_a = g_random(mean_a, dev_b) + gamma_z_a = g_random(mean_a, dev_b) + thetax(is, alive) = beta(this_box)*rott_max*Mxx(is, alive) + gamma_x_a + thetay(is, alive) = beta(this_box)*rott_max*Myy(is, alive) + gamma_y_a + thetaz(is, alive) = beta(this_box)*rott_max*Mzz(is, alive) + gamma_z_a + S_is = S_is + (beta(this_box)*rott_max*Mxx(is,alive) - thetax(is,alive))*& + (beta(this_box)*rott_max*Mxx(is,alive) - thetax(is,alive)) + S_is = S_is + (beta(this_box)*rott_max*Myy(is,alive) - thetay(is,alive))*& + (beta(this_box)*rott_max*Myy(is,alive) - thetay(is,alive)) + S_is = S_is + (beta(this_box)*rott_max*Mzz(is,alive) - thetaz(is,alive))*& + (beta(this_box)*rott_max*Mzz(is,alive) - thetaz(is,alive)) + CALL Rotate_Molecule_Axis_MPM(alive, is, thetax(is, alive), thetay(is, alive), thetaz(is, alive)) !Rigid rotate molecule around COM, the rotation angle is biased by the torque on the molecule + END DO LOOP34 + S_old_new = S_is/rott_max + S_old_new + END DO LOOP33 + + S_old_new = S_old_new /4.D0*(-1.D0) + + ! Check possible overlap, calculate the new E_vdw and E_qq, Caclulate force of new configuration, Fxx, Fyy, Fzz + + inter_overlap = .FALSE. + + LOOP35: DO is = 1, nspecies + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP SCHEDULE(STATIC) & + !$OMP PRIVATE(im, alive, inter_overlap) & + !$OMP REDUCTION(.OR.:total_overlap) + LOOP36: DO im = 1, nmolecules(is) !nmols(is,this_box) + alive = locate(im,is) + if(molecule_list(alive,is)%which_box /= this_box) cycle + if(.NOT. molecule_list(alive,is)%live) cycle + CALL Fold_Molecule(alive,is,this_box) + CALL compute_molecule_inter_overlap(alive,is,this_box,inter_overlap) + IF (inter_overlap) THEN + total_overlap = .TRUE. + END IF + END DO LOOP36 + !$OMP END PARALLEL DO + END DO LOOP35 + + IF (.NOT. total_overlap) THEN + + Conv_Drude = .TRUE. + IF(shell_mpm) CALL shell_relax(this_box, conv_drude) ! perform electrostatic energy minimization + + E_vdw_move = 0.0; E_qq_move = 0.0 + E_intra_vdw_move = 0.0; E_intra_qq_move = 0.0 + E_bond_move = 0.0 + + IF (conv_drude) THEN + + energy_cal = .TRUE. + DO is = 1, nspecies + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP SCHEDULE(STATIC) & + !$OMP PRIVATE(im, alive) & + !$OMP PRIVATE(E_vdw_move_alive, E_qq_move_alive,inter_overlap) & + !$OMP PRIVATE(E_intra_vdw_move_alive, E_intra_qq_move_alive,E_bond_alive,intra_overlap) & + !$OMP REDUCTION(+:E_vdw_move,E_qq_move, E_intra_qq_move, E_intra_vdw_move, E_bond_move) + DO im = 1, nmolecules(is) + alive = locate(im,is) + if(molecule_list(alive,is)%which_box /= this_box) cycle + if(.NOT. molecule_list(alive,is)%live) cycle + CALL Compute_Molecule_Nonbond_Inter_Force(alive, is,Fxx(is, alive, :), Fyy(is, alive, :), Fzz(is, alive, :), & + & E_vdw_move_alive, E_qq_move_alive, inter_overlap, energy_cal) + CALL Compute_Molecule_Torque(alive, is, this_box, Fxx(is, alive,:), Fyy(is, alive,:), Fzz(is, alive,:), Mxx(is, alive), Myy(is, alive), Mzz(is,alive) ) + CALL Compute_Molecule_Nonbond_Intra_Energy(alive,is,E_intra_vdw_move_alive,E_intra_qq_move_alive,intra_overlap) + CALL Compute_Molecule_Bond_Energy(alive,is,E_bond_alive) + E_vdw_move = E_vdw_move_alive + E_vdw_move + E_qq_move = E_qq_move_alive + E_qq_move + E_intra_qq_move = E_intra_qq_move_alive + E_intra_qq_move + E_intra_vdw_move = E_intra_vdw_move_alive + E_intra_vdw_move + E_bond_move = E_bond_move + E_bond_alive + END DO + !$OMP END PARALLEL DO + END DO + + + LOOP37: DO is = 1, nspecies + + S_is = 0.0 + rott_max = max_rot(is,this_box) + ! if(is .gt. 2) cycle + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP SCHEDULE(STATIC) & + !$OMP PRIVATE(im, alive) & + !$OMP REDUCTION(+:S_is) + LOOP38: DO im = 1, nmolecules(is) + alive = locate(im,is) + if(molecule_list(alive,is)%which_box /= this_box) cycle + if(.NOT. molecule_list(alive,is)%live) cycle + S_is = S_is + (beta(this_box)*rott_max*Mxx(is,alive) + thetax(is,alive))*& + (beta(this_box)*rott_max*Mxx(is,alive) + thetax(is,alive)) + S_is = S_is + (beta(this_box)*rott_max*Myy(is,alive) + thetay(is,alive))*& + (beta(this_box)*rott_max*Myy(is,alive) + thetay(is,alive)) + S_is = S_is + (beta(this_box)*rott_max*Mzz(is,alive) + thetaz(is,alive))*& + (beta(this_box)*rott_max*Mzz(is,alive) + thetaz(is,alive)) + END DO LOOP38 + !$OMP END PARALLEL DO + S_new_old = S_is/rott_max + S_new_old + END DO LOOP37 + + S_new_old = S_new_old / 4.0 * (-1.0) + + IF (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian .and. conv_drude) THEN + Energy_old_reciprocal = energy(this_box)%ewald_reciprocal + CALL Compute_System_Ewald_Reciprocal_Energy2(this_box) + END IF + + END IF + + ELSE + + DO is = 1, nspecies + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP SCHEDULE(STATIC) & + !$OMP PRIVATE(im, alive) + DO im = 1, nmolecules(is) + alive = locate(im,is) + if(molecule_list(alive,is)%which_box /= this_box) cycle + if(.NOT. molecule_list(alive,is)%live) cycle + + CALL Revert_Old_Cartesian_Coordinates(alive,is) + + END DO + !$OMP END PARALLEL DO + END DO + + END IF + + delta_e = ( E_vdw_move/2.0 - E_vdw ) + ( E_qq_move/2.0 - E_qq ) + delta_e = delta_e + E_bond_move - E_bond + E_intra_vdw_move + E_intra_qq_move - E_intra_qq - E_intra_vdw + delta_e = delta_e + energy(this_box)%ewald_reciprocal - Energy_old_reciprocal + p_acc = MIN( 1.0_DP, attempt_p * ( DEXP(-beta(this_box) * delta_e + S_new_old - S_old_new ) ) ) + + IF ( rranf() <= p_acc .and. .NOT. total_overlap .and. conv_drude ) THEN + energy(this_box)%inter_vdw = E_vdw_move/2.0 + energy(this_box)%inter_q = E_qq_move/2.0 + energy(this_box)%intra_q = E_intra_qq_move + energy(this_box)%intra_vdw = E_intra_vdw_move + energy(this_box)%bond = E_bond_move + energy(this_box)%intra = E_bond_move + energy(this_box)%total = energy(this_box)%total + delta_e + nsuccess(:,this_box)%rotation = nsuccess(:,this_box)%rotation + 1 + nequil_success(:,this_box)%rotation = nequil_success(:,this_box)%rotation + 1 + ELSE + + IF ( .NOT. total_overlap) THEN + + DO is = 1, nspecies + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP SCHEDULE(STATIC) & + !$OMP PRIVATE(im, alive) + DO im = 1, nmolecules(is) !nmols(is,this_box) + alive = locate(im,is) + if(molecule_list(alive,is)%which_box /= this_box) cycle + if(.NOT. molecule_list(alive,is)%live) cycle + CALL Revert_Old_Cartesian_Coordinates(alive,is) + END DO + !$OMP END PARALLEL DO + END DO + + IF (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian .and. conv_drude ) THEN + Energy(this_box)%ewald_reciprocal = Energy_old_reciprocal + END IF + END IF + END IF + + +IF (MOD(ntrials(1,this_box)%rotation, nupdate) == 0) THEN + DO is = 1,nspecies + IF (int_run_style == run_equil) THEN + success_ratio = REAL(nequil_success(is,this_box)%rotation,DP)/REAL(nupdate,DP) + ELSE + success_ratio = REAL(nsuccess(is,this_box)%rotation,DP)/REAL(ntrials(is,this_box)%rotation,DP) + END IF + WRITE(logunit,*) + WRITE(logunit,'(A,2x,I2,2X,A,2X,I2,2X,A,2X,F8.5)') 'Success ratio for rotation of species', is,'for box',this_box, 'is:', success_ratio + WRITE(logunit,*) 'adjust parameter B', max_rot(is,this_box) + WRITE(logunit,*) + IF (int_run_style == run_equil) THEN + nequil_success(is,this_box)%rotation = 0 + ! update the width of the species for equilibration run + IF ( success_ratio < 0.01 ) THEN + max_rot(is,this_box) = 0.1_DP*max_rot(is,this_box) + ELSE + max_rot(is,this_box) = 3.0_DP*success_ratio*max_rot(is,this_box) + END IF + END IF + ENDDO + +END IF + +CONTAINS + + SUBROUTINE Rotate_Molecule_Axis_MPM(alive, is, u, v, w) + IMPLICIT NONE + INTEGER, INTENT(IN) :: alive, is + REAL(DP), INTENT(IN) :: u, v, w + INTEGER :: axis, ia + REAL(DP) :: a, b, c + REAL(DP) :: theta, L, Lsr + REAL(DP) :: x1, y1, z1 + REAL(DP) :: x2, y2, z2 + theta = U*u + v*v + w*w + L = theta + Lsr = sqrt(L) + theta = sqrt(theta) + a = molecule_list(alive,is)%xcom + b = molecule_list(alive,is)%ycom + c = molecule_list(alive,is)%zcom + DO ia = 1, natoms(is) + x1 = atom_list(ia,alive,is)%rxp + y1 = atom_list(ia,alive,is)%ryp + z1 = atom_list(ia,alive,is)%rzp + x2 = (a*(v*v + w*w) - u*(b*v+c*w-u*x1-v*y1-w*z1))*(1.0 - cos(theta)) + L*x1*cos(theta) & + & + Lsr*(-c*v+b*w-w*y1+v*z1)*sin(theta) + y2 = (b*(u*u + w*w) - v*(a*u+c*w-u*x1-v*y1-w*z1))*(1.0 - cos(theta)) + L*y1*cos(theta) & + & + Lsr*(c*u-a*w+w*x1-u*z1)*sin(theta) + z2 = (c*(u*u + v*v) - w*(a*u+b*v-u*x1-v*y1-w*z1))*(1.0 - cos(theta)) + L*z1*cos(theta) & + & + Lsr*(-b*u+a*v-v*x1+u*y1)*sin(theta) + x2 = x2/L + y2 = y2/L + z2 = z2/L + atom_list(ia,alive,is)%rxp = x2 + atom_list(ia,alive,is)%ryp = y2 + atom_list(ia,alive,is)%rzp = z2 + END DO + END SUBROUTINE Rotate_Molecule_Axis_MPM +END SUBROUTINE Rotate_MP + + + +! Generate random number with Gaussian distribution +function g_random(m, s) +USE Type_Definitions +USE Random_Generators +implicit none +real(kind=8) :: g_random +real(kind=8) :: m , s +real(kind=8) :: x1, x2, w, y1 +real(kind=8), save :: y2 +integer(kind=4), save :: use_last = 0 + +IF(use_last == 1) THEN + y1 = y2 + use_last = 0 +ELSE + w = 2.0 + DO WHILE (w >= 1.0) + x1 = 2.0*rranf() - 1.0 + x2 = 2.0*rranf() - 1.0 + w = x1*x1 + x2*x2 + END DO + w = sqrt(-2.0*log(w)/w) + y1 = x1*w + y2 = x2*w + use_last = 1 +END IF +g_random = m + y1*s +return + +END function g_random diff --git a/Src/nptmc_driver.f90 b/Src/nptmc_driver.f90 index 93a0d645..28148204 100644 --- a/Src/nptmc_driver.f90 +++ b/Src/nptmc_driver.f90 @@ -59,7 +59,7 @@ SUBROUTINE NPTMC_Driver IMPLICIT NONE -! !$ include 'omp_lib.h' +! !$ include 'omp_lib.h' INTEGER :: i,j,k, this_box, ibox, is, ifrag, which_step, ireac INTEGER :: howmanyfrac, ii,jj, im1, im2, alive1, alive2 @@ -135,7 +135,6 @@ SUBROUTINE NPTMC_Driver i = i + 1 ! We will select a move from Golden Sampling scheme - rand_no = rranf() which_step = i @@ -146,8 +145,12 @@ SUBROUTINE NPTMC_Driver ELSE !$ time_s = omp_get_wtime() END IF - - CALL Translate(this_box,which_step) + + IF(shell_mpm) THEN + CALL Translate_MP(this_box,which_step) + ELSE + CALL Translate(this_box,which_step) + END IF IF(.NOT. openmp_flag) THEN CALL cpu_time(time_e) @@ -164,8 +167,12 @@ SUBROUTINE NPTMC_Driver ELSE !$ time_s = omp_get_wtime() END IF - - CALL Rotate(this_box) + + IF(shell_mpm) THEN + CALL Rotate_MP(this_box,which_step) + ELSE + CALL Rotate(this_box) + END IF IF(.NOT. openmp_flag) THEN CALL cpu_time(time_e) diff --git a/Src/nvtmc_driver.f90 b/Src/nvtmc_driver.f90 index b12ae811..40a5304c 100644 --- a/Src/nvtmc_driver.f90 +++ b/Src/nvtmc_driver.f90 @@ -86,7 +86,7 @@ SUBROUTINE NVTMC_Driver openmp_flag = .FALSE. i = 0 - + this_box = 1 !$ openmp_flag = .TRUE. write(*,*) 'openmp_flag = ', openmp_flag @@ -104,8 +104,8 @@ SUBROUTINE NVTMC_Driver ! We will select a move from Golden Sampling scheme - rand_no = rranf() + IF (rand_no <= cut_trans) THEN @@ -115,7 +115,13 @@ SUBROUTINE NVTMC_Driver !$ time_s = omp_get_wtime() END IF - CALL Translate(this_box,which_step) + IF(shell_mpm) THEN + CALL Translate_MP(this_box,which_step) + ELSE + CALL Translate(this_box,which_step) + END IF + + !CALL Translate(this_box,which_step) IF(.NOT. openmp_flag) THEN CALL cpu_time(time_e) @@ -132,8 +138,14 @@ SUBROUTINE NVTMC_Driver ELSE !$ time_s = omp_get_wtime() END IF - - CALL Rotate(this_box) + + IF(shell_mpm) THEN + CALL Rotate_MP(this_box,which_step) + ELSE + CALL Rotate(this_box) + END IF + + !CALL Rotate(this_box) IF(.NOT. openmp_flag) THEN CALL cpu_time(time_e) diff --git a/Src/precalculate.f90 b/Src/precalculate.f90 index b79e829c..ab325d44 100644 --- a/Src/precalculate.f90 +++ b/Src/precalculate.f90 @@ -100,3 +100,109 @@ SUBROUTINE precalculate ALLOCATE(energy(nbr_boxes),virial(nbr_boxes)) END SUBROUTINE precalculate + + + + +SUBROUTINE generate_reaf_table(ibox) + +USE Run_variables +USE Energy_Routines + +IMPLICIT NONE + +integer, intent(in) :: ibox +integer :: i, j, k, bin +real :: cut_reaf +real :: r_rc, p_rc, p_rc_prime, du_dr, rij, rijsq +cut_reaf = rcut_coul(ibox) + +DO i = 1, int(cut_reaf/reaf_width) + 1 + +rij = i*reaf_width +rijsq = rij*rij +r_rc = rij/rcut_coul(ibox) + +p_rc = (1.0 - r_rc) +p_rc_prime = -4.0*p_rc*p_rc*p_rc*(1.0+1.6*r_rc+0.4*r_rc*r_rc) +p_rc_prime = p_rc_prime + p_rc*p_rc*p_rc*p_rc*(1.6+0.8*r_rc) +p_rc = p_rc*p_rc*p_rc*p_rc*(1.0 + 1.6_DP*r_rc + 0.4_DP*r_rc*r_rc) +dU_dr = (-1.0)/rijsq*p_rc + 1.0/sqrt(rijsq*rcut_coulsq(ibox))*p_rc_prime +du_dr = du_dr*(-1.0)*charge_factor + + +reaf(i,ibox) = dU_dr + +END DO + +END SUBROUTINE generate_reaf_table + +SUBROUTINE Initialize_MP + +USE Type_Definitions +USE Run_variables +USE Energy_Routines + +!LOGICAL, INTENT(OUT), optional :: mpm_logical + + shell_mpm = .FALSE. + + DO i = 1, nbr_boxes + DO is = 1, nspecies + DO im = 1, nmolecules(is) + DO ia = 1, natoms(is) + atom_list(ia,im,is)%drude_type = .FALSE. + if(nonbond_list(ia,is)%element == 'G' ) then + atom_list(ia,im,is)%drude_type = .TRUE. + shell_mpm = .TRUE. + end if + END DO + END DO + END DO + END DO + + ! mpm force and torque + + IF(shell_mpm) THEN + + ALLOCATE( Mxx(nspecies, MAXVAL(nmolecules)), stat = AllocateStatus) + ALLOCATE( Myy(nspecies, MAXVAL(nmolecules)), stat = AllocateStatus) + ALLOCATE( Mzz(nspecies, MAXVAL(nmolecules)), stat = AllocateStatus) + + ALLOCATE( Fxx(nspecies, MAXVAL(nmolecules), MAXVAL(natoms)) , stat = AllocateStatus) + ALLOCATE( Fyy(nspecies, MAXVAL(nmolecules), MAXVAL(natoms)) , stat = AllocateStatus) + ALLOCATE( Fzz(nspecies, MAXVAL(nmolecules), MAXVAL(natoms)) , stat = AllocateStatus) + + ALLOCATE(drudeFx(nspecies,maxval(nmolecules),maxval(natoms))) + ALLOCATE(drudeFy(nspecies,maxval(nmolecules),maxval(natoms))) + ALLOCATE(drudeFz(nspecies,maxval(nmolecules),maxval(natoms))) + ALLOCATE(dmpmx(nspecies,SUM(nmolecules)), dmpmy(nspecies,sum(nmolecules)), dmpmz(nspecies,sum(nmolecules)) ) + ALLOCATE(thetax(nspecies,SUM(nmolecules)), thetay(nspecies,sum(nmolecules)), thetaz(nspecies,sum(nmolecules))) + + drudeFx = 0.0_DP + drudeFy = 0.0_DP + drudeFz = 0.0_DP + + Fxx = 0.0; Fyy = 0.0 ; Fzz = 0.0 + + Mxx = 0.0_DP; Myy = 0.0_DP; Mzz = 0.0_DP + + dmpmx = 0.0_DP; dmpmy = 0.0_DP; dmpmz = 0.0_DP + thetax = 0.0_DP; thetay = 0.0_DP; thetaz = 0.0_DP + + reaf_width = 0.01 + + ALLOCATE(reaf(int(maxval(rcut_coul)/reaf_width)+1,nbr_boxes)) + + reaf = 0.0_DP + + ! Generate table for reaction field force + DO ibox = 1, nbr_boxes + CALL generate_reaf_table(ibox) + END DO + + mpm_logical = .TRUE. + +END IF + +END SUBROUTINE Initialize_MP diff --git a/Src/run_variables.f90 b/Src/run_variables.f90 index 5a047de8..0755d996 100644 --- a/Src/run_variables.f90 +++ b/Src/run_variables.f90 @@ -104,12 +104,15 @@ MODULE Run_Variables INTEGER, PARAMETER :: vdw_charmm = 6 INTEGER, PARAMETER :: vdw_cut_switch = 7 INTEGER, PARAMETER :: vdw_mie = 8 + INTEGER, PARAMETER :: born_cut_tail = 9 INTEGER, PARAMETER :: charge_none = 0 INTEGER, PARAMETER :: charge_coul = 1 INTEGER, PARAMETER :: charge_cut = 2 INTEGER, PARAMETER :: charge_ewald = 3 INTEGER, PARAMETER :: charge_minimum = 4 + INTEGER, PARAMETER :: charge_gaussian = 5 + REAL(DP), DIMENSION(:), ALLOCATABLE :: rcut_cbmc REAL(DP), DIMENSION(:), ALLOCATABLE :: rcut_vdw, rcut_coul, ron_charmm, roff_charmm, rcut_max @@ -126,7 +129,7 @@ MODULE Run_Variables INTEGER, DIMENSION(:,:), ALLOCATABLE ::vdw_int_table REAL(DP), DIMENSION(:,:), ALLOCATABLE :: vdw_param1_table REAL(DP), DIMENSION(:,:), ALLOCATABLE :: vdw_param2_table, vdw_param3_table - REAL(DP), DIMENSION(:,:), ALLOCATABLE :: vdw_param4_table, vdw_param5_table + REAL(DP), DIMENSION(:,:), ALLOCATABLE :: vdw_param4_table, vdw_param5_table!, vdw_param6_table REAL(DP), DIMENSION(:), ALLOCATABLE :: alpha_ewald, h_ewald_cut REAL(DP), DIMENSION(:), ALLOCATABLE :: alphal_ewald REAL(DP), DIMENSION(:), ALLOCATABLE :: ewald_p_sqrt, ewald_p @@ -518,6 +521,10 @@ MODULE Run_Variables REAL(DP), DIMENSION(:), ALLOCATABLE :: P_inst REAL(DP), DIMENSION(:), ALLOCATABLE :: P_ideal + REAL(DP), DIMENSION(:), ALLOCATABLE :: P_inst_plus, P_inst_minus + + LOGICAL :: constant_vol, P_calc, p_plus_1 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Energy check !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -544,8 +551,32 @@ MODULE Run_Variables !!!! Zeolite variables -REAL(DP), ALLOCATABLE, DIMENSION(:) :: x_lat, y_lat, z_lat -INTEGER :: n_lat_atoms - + REAL(DP), ALLOCATABLE, DIMENSION(:) :: x_lat, y_lat, z_lat + INTEGER :: n_lat_atoms + +!Gaussian charge variables + REAL(DP), ALLOCATABLE :: alp_ij(:,:,:,:), alp_ia(:,:,:) +! Reaction field table + REAL(DP), DIMENSION(:,:) , ALLOCATABLE :: reaf + REAL(DP) :: reaf_width + +! shell particle + LOGICAL :: shell_mpm + +! multiparticle move + REAL(DP), DIMENSION(:,:,:), ALLOCATABLE :: Fxx, Fyy, Fzz + REAL(DP), DIMENSION(:,:,:), ALLOCATABLE :: Fxx_old, Fyy_old, Fzz_old + REAL(DP), DIMENSION(:,:), ALLOCATABLE :: Mxx, Myy, Mzz + REAL(DP), DIMENSION(:,:), ALLOCATABLE :: Mxx_old, Myy_old, Mzz_old + REAL(DP), DIMENSION(:,:,:), ALLOCATABLE :: drudeFx, drudeFy, drudeFz + + REAL(DP), DIMENSION(:,:), ALLOCATABLE :: dmpmx,dmpmy,dmpmz + REAL(DP), DIMENSION(:,:), ALLOCATABLE :: thetax,thetay,thetaz + + INTEGER :: check_err ! counter for unconverged EM + + ! pressure calculation using trial volume change in GEMC-NVT simulation + + END MODULE Run_Variables diff --git a/Src/shell_gemc_nvt_volume.f90 b/Src/shell_gemc_nvt_volume.f90 new file mode 100644 index 00000000..921b8638 --- /dev/null +++ b/Src/shell_gemc_nvt_volume.f90 @@ -0,0 +1,688 @@ +!******************************************************************************** +! Cassandra - An open source atomistic Monte Carlo software package +! developed at the University of Notre Dame. +! http://cassandra.nd.edu +! Prof. Edward Maginn +! Copyright (2013) University of Notre Dame du Lac +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see . +!******************************************************************************* + +SUBROUTINE Shell_GEMC_NVT_Volume(box1, box2) + + !*********************************************************************** + ! + ! The routine performs an NVT volume move for a GEMC NVT simulation. There are core-shell units in the box. + ! At present the simulation is restricted to two boxes. Future possibility + ! may include inclusion of multiple boxes. + ! + ! Called by + ! + ! gemc_driver + ! + ! ! + !************************************************************************* + + USE Run_Variables + USE Random_Generators + USE Volume + USE IO_Utilities + USE Energy_Routines + USE Pair_Nrg_Routines + + IMPLICIT NONE + + + INTEGER, INTENT(OUT) :: box1, box2 + + INTEGER :: tot_mol_box1, tot_mol_box2, first_box, second_box + INTEGER :: ibox, nvecs_old_1, nvecs_old_2, nvecs_max + INTEGER :: is, im + + REAL(DP) :: delta_volume, factor, delta_e_1, delta_e_2 + REAL(DP) :: success_ratio + REAL(DP), DIMENSION(maxk) :: hx_old_1, hy_old_1, hz_old_1, Cn_old_1 + REAL(DP), DIMENSION(maxk) :: hx_old_2, hy_old_2, hz_old_2, Cn_old_2 + REAL(DP) :: v_ratio_o, v_total, vol_factor + + LOGICAL :: overlap, accept, accept_or_reject, allocation_cos_sin + + TYPE(Box_Class) :: box_list_old_1, box_list_old_2 + TYPE(Energy_Class) :: energy_old_1, energy_old_2 + + REAL(DP) :: checke + LOGICAL :: superbad + + INTEGER :: position, alive, my_box + + REAL(DP) :: rcut_vdw_old_1, rcut_coul_old_1, rcut3_old_1, rcut9_old_1, alpha_ewald_old_1 + REAL(DP) :: h_ewald_cut_old_1, rcut_vdwsq_old_1, rcut_coulsq_old_1, rcut_vdw3_old_1 + REAL(DP) :: rcut_vdw6_old_1, rcut_max_old_1 + + REAL(DP) :: rcut_vdw_old_2, rcut_coul_old_2, rcut3_old_2, rcut9_old_2, alpha_ewald_old_2 + REAL(DP) :: h_ewald_cut_old_2, rcut_vdwsq_old_2, rcut_coulsq_old_2, rcut_vdw3_old_2 + REAL(DP) :: rcut_vdw6_old_2, rcut_max_old_2 + REAL(DP) :: E_vdw_move_alive, E_qq_move_alive + INTEGER :: i, j, box_number + LOGICAL :: conv_drude1, conv_drude2 + + IF (f_dv) THEN + ! Pick a box at random + box1 = INT (rranf() * nbr_boxes) + 1 + + ! Obtain the identity of the other box + + DO WHILE(.TRUE.) + + box2 = INT ( rranf() * nbr_boxes) + 1 + + IF ( box2 /= box1 ) EXIT + + END DO + + ELSE IF (f_vratio) THEN + + box1 = 1 + box2 = 2 + + END IF + + IF(constant_vol) THEN !perform trial volume change to calculate pressure + p_calc = .FALSE. + IF(P_plus_1) THEN + box1 = 1 + box2 = 2 + ELSE + box1 = 2 + box2 = 1 + END IF + END IF + + IF(.NOT. constant_vol) THEN + + tot_trials(box1) = tot_trials(box1) + 1 + tot_trials(box2) = tot_trials(box2) + 1 + + nvolumes(box1) = nvolumes(box1) + 1 + nvolumes(box2) = nvolumes(box2) + 1 + END IF + + ! store old cell matrix + + box_list_old_1 = box_list(box1) + box_list_old_2 = box_list(box2) + + ! Store the old configurations of all atoms and COMs + + CALL Save_Cartesian_Coordinates_Box(box1,tot_mol_box1) + CALL Save_Cartesian_Coordinates_Box(box2,tot_mol_box2) + + ! DO NOT store the pair interactions + + ! this performs a trial volume change move to calculate pressure (without evalulation of virial) + ! Efficient pressure estimation in molecular simulations without evaluating the virial. J. Che. Phys. 105, 8469 + ! P = kb*T/delta_Volume * ln< (V'/V)^N * exp(-beta*delta_Energy) > + ! The output is the < (V'/V)^N * exp(-beta*delta_Energy) > + ! The system pressure needs to be calculated outside Cassandra + + IF (f_dv) THEN + + IF(constant_vol) THEN + delta_volume = 25.0 !box_list(box1)%dv_max + ! print *, delta_volume + ELSE + delta_volume = (1.0_DP - 2.0_DP * rranf()) * box_list(box1)%dv_max + END IF + + ! update each of the box volumes + + box_list(box1)%volume = box_list(box1)%volume + delta_volume + box_list(box2)%volume = box_list(box2)%volume - delta_volume + + ELSE IF (f_vratio) THEN + + print *, 'ratio' + v_ratio_o = (box_list(box1)%volume/box_list(box2)%volume) + v_total = box_list(box1)%volume + box_list(box2)%volume + + delta_volume = (1.0_DP - 2.0_DP * rranf()) * box_list(box1)%dv_max + vol_factor = v_ratio_o * EXP(delta_volume) + box_list(box1)%volume = v_total * vol_factor / (1.0_DP + vol_factor) + box_list(box2)%volume = v_total - box_list(box1)%volume + + END IF + + IF ( box_list(box1)%int_box_shape == int_cubic ) THEN + + box_list(box1)%length(1,1) = (box_list(box1)%volume)**(1.0_DP/3.0_DP) + box_list(box1)%length(2,2) = box_list(box1)%length(1,1) + box_list(box1)%length(3,3) = box_list(box1)%length(2,2) + + ELSE + + err_msg = '' + err_msg = 'Noncubic shape encountered for the box' + err_msg = Int_To_String(box1) + CALL Clean_Abort(err_msg, 'GEMC_NVT_Volume') + + END IF + + IF ( box_list(box2)%int_box_shape == int_cubic ) THEN + + box_list(box2)%length(1,1) = (box_list(box2)%volume) ** (1.0_DP/3.0_DP) + box_list(box2)%length(2,2) = box_list(box2)%length(1,1) + box_list(box2)%length(3,3) = box_list(box2)%length(2,2) + + ELSE + + err_msg = '' + err_msg = 'Noncubic shape encountered for the box' + err_msg = Int_To_String(box2) + CALL Clean_Abort(err_msg, 'GEMC_NVT_Volume') + + END IF + + CALL Compute_Cell_Dimensions(box1) + CALL Compute_Cell_Dimensions(box2) + + IF ( l_half_len_cutoff(box1)) THEN + IF ( box_list(box1)%int_box_shape == int_cubic ) THEN + + ! store old cutoffs and other associated quantities + + rcut_vdw_old_1 = rcut_vdw(box1) + rcut_coul_old_1 = rcut_coul(box1) + rcut_vdwsq_old_1 = rcut_vdwsq(box1) + rcut_coulsq_old_1 = rcut_coulsq(box1) + + rcut3_old_1 = rcut3(box1) + rcut9_old_1 = rcut9(box1) + rcut_vdw3_old_1 = rcut_vdw3(box1) + rcut_vdw6_old_1 = rcut_vdw6(box1) + + rcut_max_old_1 = rcut_max(box1) + + IF (int_charge_sum_style(box1) == charge_ewald .or. int_charge_sum_style(box1) == charge_gaussian) THEN + +! alpha_ewald_old_1 = alpha_ewald(box1) + h_ewald_cut_old_1 = h_ewald_cut(box1) + + END IF + + rcut_vdw(box1) = 0.5_DP * box_list(box1)%length(1,1) + rcut_coul(box1) = rcut_vdw(box1) + rcut_vdwsq(box1) = rcut_vdw(box1) * rcut_vdw(box1) + rcut_coulsq(box1) = rcut_vdwsq(box1) + + rcut_vdw3(box1) = rcut_vdwsq(box1) * rcut_vdw(box1) + rcut_vdw6(box1) = rcut_vdw3(box1) * rcut_vdw3(box1) + rcut3(box1) = rcut_vdw3(box1) + rcut9(box1) = rcut3(box1) * rcut_vdw6(box1) + + rcut_max(box1) = rcut_vdw(box1) + + IF (int_charge_sum_style(box1) == charge_ewald .or. int_charge_sum_style(box1) == charge_gaussian) THEN + +! alpha_ewald(box1) = ewald_p_sqrt(box1) / rcut_coul(box1) + h_ewald_cut(box1) = 2.0_DP * ewald_p(box1) / rcut_coul(box1) + + END IF + + END IF + + ELSE + + IF ( box_list(box1)%int_box_shape == int_cubic ) THEN + IF ( 0.5_DP * box_list(box1)%length(1,1) < rcut_vdw(box1) .OR. & + 0.5_DP * box_list(box1)%length(1,1) < rcut_coul(box1) .OR. & + 0.5_DP * box_list(box1)%length(1,1) < roff_charmm(box1) .OR. & + 0.5_DP * box_list(box1)%length(1,1) < roff_switch(box1) ) THEN + err_msg = '' + err_msg(1) = 'Cutoff is greater than the half box length' + err_msg(2) = Int_to_String(box1) + CALL Clean_Abort(err_msg,'GEMC_NVT_Volume') + END IF + END IF + + END IF + + IF ( l_half_len_cutoff(box2)) THEN + IF ( box_list(box2)%int_box_shape == int_cubic ) THEN + + ! store old cutoffs and other associated quantities + + rcut_vdw_old_2 = rcut_vdw(box2) + rcut_coul_old_2 = rcut_coul(box2) + rcut_vdwsq_old_2 = rcut_vdwsq(box2) + rcut_coulsq_old_2 = rcut_coulsq(box2) + + rcut3_old_2 = rcut3(box2) + rcut9_old_2 = rcut9(box2) + rcut_vdw3_old_2 = rcut_vdw3(box2) + rcut_vdw6_old_2 = rcut_vdw6(box2) + + rcut_max_old_2 = rcut_max(box2) + + IF (int_charge_sum_style(box2) == charge_ewald .or. int_charge_sum_style(box2) == charge_gaussian) THEN + +! alpha_ewald_old_2 = alpha_ewald(box2) + h_ewald_cut_old_2 = h_ewald_cut(box2) + + END IF + + rcut_vdw(box2) = 0.5_DP * box_list(box2)%length(1,1) + rcut_coul(box2) = rcut_vdw(box2) + rcut_vdwsq(box2) = rcut_vdw(box2) * rcut_vdw(box2) + rcut_coulsq(box2) = rcut_vdwsq(box2) + + rcut_vdw3(box2) = rcut_vdwsq(box2) * rcut_vdw(box2) + rcut_vdw6(box2) = rcut_vdw3(box2) * rcut_vdw3(box2) + rcut3(box2) = rcut_vdw3(box2) + rcut9(box2) = rcut3(box2) * rcut_vdw6(box2) + + rcut_max(box2) = rcut_vdw(box2) + + IF (int_charge_sum_style(box2) == charge_ewald .or. int_charge_sum_style(box2) == charge_gaussian) THEN + +! alpha_ewald(box2) = ewald_p_sqrt(box2) / rcut_coul(box2) + h_ewald_cut(box2) = 2.0_DP * ewald_p(box2) / rcut_coul(box2) + + END IF + + END IF + + + ELSE + + IF ( box_list(box2)%int_box_shape == int_cubic ) THEN + IF ( 0.5_DP * box_list(box2)%length(1,1) < rcut_vdw(box2) .OR. & + 0.5_DP * box_list(box2)%length(1,1) < rcut_coul(box2) .OR. & + 0.5_DP * box_list(box2)%length(1,1) < roff_charmm(box2) .OR. & + 0.5_DP * box_list(box2)%length(1,1) < roff_switch(box2) ) THEN + err_msg = '' + err_msg(1) = 'Cutoff is greater than the half box length' + err_msg(2) = Int_To_String(box2) + CALL Clean_Abort(err_msg,'Shell_GEMC_NVT_Volume') + END IF + END IF + + END IF + + + ! Rescale the COM and all atomic positions + + CALL Scale_COM_Cartesian(box1,box_list_old_1) + CALL Scale_COM_Cartesian(box2,box_list_old_2) + + ! Now let us compute the energy change due to the combined move + + energy_old_1 = energy(box1) + energy_old_2 = energy(box2) + + IF ( int_charge_sum_style(box1) == charge_ewald .or. int_charge_sum_style(box1) == charge_gaussian) THEN + + nvecs_old_1 = nvecs(box1) + nvecs_old_2 = nvecs(box2) + nvecs_max = MAXVAL(nvecs) + + !!$OMP PARALLEL WORKSHARE DEFAULT(SHARED) + + hx_old_1(:) = hx(:,box1) + hy_old_1(:) = hy(:,box1) + hz_old_1(:) = hz(:,box1) + Cn_old_1(:) = Cn(:,box1) + + hx_old_2(:) = hx(:,box2) + hy_old_2(:) = hy(:,box2) + hz_old_2(:) = hz(:,box2) + Cn_old_2(:) = Cn(:,box2) + + !!$OMP END PARALLEL WORKSHARE + + ! Determine the new k vectors for this box. The call will change Cn, hx, hy and hz and hence will + ! change cos_sum and sin_sum. + + CALL Ewald_Reciprocal_Lattice_Vector_Setup(box1) + CALL Ewald_Reciprocal_Lattice_Vector_Setup(box2) + + ! check to see if we need to reallocate arrays + + allocation_cos_sin = .FALSE. + + IF ((nvecs(box1) > SIZE(cos_sum,1)) .OR. (nvecs(box2) > SIZE(cos_sum,1))) THEN + ! we need to expand the cos_sum and sin_sum arrays + allocation_cos_sin = .TRUE. + + DEALLOCATE(cos_sum,sin_sum) + + IF (ALLOCATED(cos_mol)) DEALLOCATE(cos_mol) + IF (ALLOCATED(sin_mol)) DEALLOCATE(sin_mol) + ALLOCATE(cos_sum(MAXVAL(nvecs),nbr_boxes),Stat=AllocateStatus) + + IF (AllocateStatus /= 0) THEN + err_msg = '' + err_msg(1) = 'Memory could not be allocated for cos_sum' + err_msg(2) = 'allocation_cos_sin' + CALL Clean_Abort(err_msg,'GEMC_NVT_VOLUME') + END IF + + + ALLOCATE(sin_sum(MAXVAL(nvecs),nbr_boxes), Stat = AllocateStatus) + + IF (Allocatestatus /= 0) THEN + err_msg = '' + err_msg(1) = 'Memory could not be allocated for sin_sum' + err_msg(2) = 'allocation_cos_sin' + CALL Clean_Abort(err_msg,'GEMC_NVT_VOLUME') + END IF + + ALLOCATE(cos_mol(MAXVAL(nvecs),SUM(nmolecules)), Stat = AllocateStatus) + + IF (AllocateStatus /= 0) THEN + err_msg = '' + err_msg(1) = 'Memory could not be allocated for cos_mol' + err_msg(2) = 'allocation_cos_sin' + CALL Clean_Abort(err_msg,'gemc_volume_change.f90') + END IF + + ALLOCATE(sin_mol(MAXVAL(nvecs),SUM(nmolecules)), Stat = AllocateStatus) + + IF (AllocateStatus /= 0) THEN + err_msg = '' + err_msg(1) = 'Memory could not be allocated for sin_mol' + err_msg(2) = 'allocation_cos_sin' + CALL Clean_Abort(err_msg,'gemc_volume_change.f90') + END IF + + END IF + + END IF + + ! Compute the energies of the boxes. For computational efficiency, + ! we will determine the energy of the box whose volume decreased + ! as this box will likely to have more overlaps + + IF (box_list(box1)%volume < box_list_old_1%volume) THEN + ! box1 is shrinking + first_box = box1 + second_box = box2 + ELSE + ! else box2 is shrinking + first_box = box2 + second_box = box1 + END IF + + + conv_drude1 = .TRUE. + IF (shell_mpm) CALL shell_relax(first_box, conv_drude1) !perform energy minimization + IF (conv_drude1) CALL Compute_Total_System_Energy(first_box, .TRUE., overlap) + + IF (overlap .or. .NOT. conv_drude1 ) THEN + + CALL Reset_Coords + + ELSE + + conv_drude2 = .TRUE. + IF(shell_mpm) CALL shell_relax(second_box, conv_drude2) !perform energy minimization + IF (conv_drude2) CALL Compute_Total_System_Energy(second_box, .TRUE.,overlap) + + ! actually there should be no overlap for the box whose dimensions + ! increase but we will include this check only for safety. + + IF (overlap .or. .NOT. conv_drude2) THEN + CALL Reset_Coords + ELSE + ! accept or reject the move based on the acceptance rule + + delta_e_1 = energy(box1)%total - energy_old_1%total + delta_e_2 = energy(box2)%total - energy_old_2%total + + IF (f_dv) THEN + factor = beta(box1) * delta_e_1 + beta(box2) * delta_e_2 - & + tot_mol_box1 * DLOG( box_list(box1)%volume / box_list_old_1%volume) - & + tot_mol_box2 * DLOG( box_list(box2)%volume / box_list_old_2%volume) + + ELSE IF(f_vratio) THEN + + factor = beta(box1) * delta_e_1 + beta(box2) * delta_e_2 - & + REAL((tot_mol_box1 + 1),DP) * DLOG( box_list(box1)%volume / box_list_old_1%volume) - & + REAL((tot_mol_box2 + 1),DP) * DLOG( box_list(box2)%volume / box_list_old_2%volume) + + END IF + + IF (box_list(box1)%volume > box_list_old_1%volume .and. constant_vol) THEN + !print *, 'delta_volume', delta_volume, box1, box2, p_inst_plus(:) + P_inst_plus(box1) = exp(sum(nmols(:,box1)) * log(box_list(box1)%volume/box_list_old_1%volume) - beta(box1)*delta_e_1) !pressure change due to a positive volume change + P_inst_minus(box2) = exp(sum(nmols(:,box2)) * log(box_list(box2)%volume/box_list_old_2%volume) - beta(box2)*delta_e_2) !pressure change due to a negative volume change + P_calc = .TRUE. + ELSE IF(constant_vol) THEN + P_inst_minus(box1) = exp(sum(nmols(:,box1)) * log(box_list(box1)%volume/box_list_old_1%volume) - beta(box1)*delta_e_1) + P_inst_plus(box2) = exp(sum(nmols(:,box2)) * log(box_list(box2)%volume/box_list_old_2%volume) - beta(box2)*delta_e_2) + P_calc = .TRUE. + END IF + + accept = accept_or_reject(factor) + + IF (accept .and. .NOT. constant_vol) THEN + + nvol_success(box1) = nvol_success(box1) + 1 + nvol_success(box2) = nvol_success(box2) + 1 + ivol_success(box1) = ivol_success(box1) + 1 + ivol_success(box2) = ivol_success(box2) + 1 + + ! energy,positions and box dimensions are already updated + IF (int_charge_sum_style(box1) == charge_ewald .or. int_charge_sum_style(box1) == charge_gaussian) THEN + + IF (allocation_cos_sin) THEN + ! Now deallocate cos_sum_old and sin_sum_old so that they have the same dimensions + ! as sin_sum and cos_sum + DEALLOCATE(cos_sum_old,sin_sum_old) + ALLOCATE(cos_sum_old(SIZE(cos_sum,1),nbr_boxes)) + ALLOCATE(sin_sum_old(SIZE(sin_sum,1),nbr_boxes)) + END IF + END IF + + ELSE + + CALL Reset_Coords + + END IF + + END IF + + END IF + + ! Update the maximum volume modulus of equilibration runs + + IF (MOD(nvolumes(box1),nvol_update) == 0 ) THEN + + ! IF ( int_run_style == run_equil) THEN + + success_ratio = REAL(ivol_success(box1),DP)/REAL(nvol_update,DP) + ivol_success(box1) = 0 + ivol_success(box2) = 0 + + ! Do not update dv_max for the evaluation of pressure + ! A uniform dv is used for pressure evaluation + + IF (box_list(box1)%int_box_shape == int_cubic) THEN + IF ( success_ratio < 0.0001 ) THEN + ! decrease the maximum displacement by 5% + + box_list(:)%dv_max = 0.1_DP * box_list(:)%dv_max + + ELSE + + box_list(:)%dv_max = 3.0_DP*success_ratio * box_list(:)%dv_max + + END IF + + ! END IF + + ELSE + + success_ratio = REAL(nvol_success(box1),DP)/REAL(nvolumes(box1),DP) + + END IF + + WRITE(logunit,*) + WRITE(logunit,'(A35,2X,F10.4,2X,A7,2X,I2)')'Successful volume attempt ratio is ', & + success_ratio, 'for box', box1 + + END IF + + CONTAINS + + SUBROUTINE Reset_Coords + + IMPLICIT NONE + + CALL Reset_Cartesian_Coordinates_Box(box1) + CALL Reset_Cartesian_Coordinates_Box(box2) + + ! box list and energy + + box_list(box1) = box_list_old_1 + box_list(box2) = box_list_old_2 + + energy(box1) = energy_old_1 + energy(box2) = energy_old_2 + + + IF (l_half_len_cutoff(box1)) THEN + + rcut_vdw(box1) = rcut_vdw_old_1 + rcut_coul(box1) = rcut_coul_old_1 + rcut_vdwsq(box1) = rcut_vdwsq_old_1 + rcut_coulsq(box1) = rcut_coulsq_old_1 + + rcut3(box1) = rcut3_old_1 + rcut9(box1) = rcut9_old_1 + rcut_vdw3(box1) = rcut_vdw3_old_1 + rcut_vdw6(box1) = rcut_vdw6_old_1 + + rcut_max(box1) = rcut_max_old_1 + + IF( int_charge_sum_style(box1) == charge_ewald .or. int_charge_sum_style(box1) == charge_gaussian ) THEN + + !alpha_ewald(box1) = alpha_ewald_old_1 + h_ewald_cut(box1) = h_ewald_cut_old_1 + + END IF + + + END IF + + IF (l_half_len_cutoff(box2)) THEN + + rcut_vdw(box2) = rcut_vdw_old_2 + rcut_coul(box2) = rcut_coul_old_2 + rcut_vdwsq(box2) = rcut_vdwsq_old_2 + rcut_coulsq(box2) = rcut_coulsq_old_2 + + rcut3(box2) = rcut3_old_2 + rcut9(box2) = rcut9_old_2 + rcut_vdw3(box2) = rcut_vdw3_old_2 + rcut_vdw6(box2) = rcut_vdw6_old_2 + + rcut_max(box2) = rcut_max_old_2 + + IF( int_charge_sum_style(box2) == charge_ewald .or. int_charge_sum_style(box2) == charge_gaussian ) THEN + +! alpha_ewald(box2) = alpha_ewald_old_2 + h_ewald_cut(box2) = h_ewald_cut_old_2 + + END IF + + + END IF + + IF (int_charge_sum_style(box1) == charge_ewald .or. int_charge_sum_style(box1) == charge_gaussian) THEN + + nvecs(box1) = nvecs_old_1 + nvecs(box2) = nvecs_old_2 + + IF ( allocation_cos_sin) THEN + + DEALLOCATE(cos_sum,sin_sum) + DEALLOCATE(cos_mol,sin_mol) + + ALLOCATE(cos_sum(MAXVAL(nvecs),nbr_boxes),stat = AllocateStatus) + + IF (Allocatestatus /= 0) THEN + err_msg = '' + err_msg(1) = 'Memory could not be allocated for cos_sum' + err_msg(2) = 'volume move rejected' + CALL Clean_Abort(err_msg,'Volume_Change') + END IF + + ALLOCATE(sin_sum(MAXVAL(nvecs),nbr_boxes),Stat = Allocatestatus) + + IF (Allocatestatus /= 0) THEN + err_msg = '' + err_msg(1) = 'Memory could not be allocated in the volume rejection' + CALL Clean_Abort(err_msg,'Volume_Change') + END IF + + ALLOCATE(cos_mol(MAXVAL(nvecs),SUM(nmolecules)),Stat = Allocatestatus) + + IF (Allocatestatus /= 0) THEN + err_msg = '' + err_msg(1) = 'Memory could not be allocated for cos_mol in the volume rejection' + CALL Clean_Abort(err_msg,'GEMC NVT Volume_Change') + END IF + + ALLOCATE(sin_mol(MAXVAL(nvecs),SUM(nmolecules)),Stat = Allocatestatus) + + IF (Allocatestatus /= 0) THEN + err_msg = '' + err_msg(1) = 'Memory could not be allocated for sin_mol in the volume rejection' + CALL Clean_Abort(err_msg,'GEMC NVT Volume_Change') + END IF + + END IF + + !!$OMP PARALLEL WORKSHARE DEFAULT(SHARED) + + cos_sum(:,:) = cos_sum_old(:,:) + sin_sum(:,:) = sin_sum_old(:,:) + + hx(:,box1) = hx_old_1(:) + hy(:,box1) = hy_old_1(:) + hz(:,box1) = hz_old_1(:) + Cn(:,box1) = Cn_old_1(:) + + hx(:,box2) = hx_old_2(:) + hy(:,box2) = hy_old_2(:) + hz(:,box2) = hz_old_2(:) + Cn(:,box2) = Cn_old_2(:) + + !!$OMP END PARALLEL WORKSHARE + + ! DEALLOCATE(cos_mol_old,sin_mol_old) + + ! here we make sure that cos_sum_old and sin_sum_old have the same dimensions + ! as cos_sum and sin_sum + + DEALLOCATE(cos_sum_old,sin_sum_old) + ALLOCATE(cos_sum_old(SIZE(cos_sum,1),nbr_boxes),sin_sum_old(SIZE(sin_sum,1),nbr_boxes)) + + END IF + + END SUBROUTINE Reset_Coords + +END SUBROUTINE Shell_GEMC_NVT_Volume + + diff --git a/Src/shell_gemc_particle_transfer.f90 b/Src/shell_gemc_particle_transfer.f90 new file mode 100644 index 00000000..c277dd46 --- /dev/null +++ b/Src/shell_gemc_particle_transfer.f90 @@ -0,0 +1,431 @@ +!******************************************************************************** +! Cassandra - An open source atomistic Monte Carlo software package +! developed at the University of Notre Dame. +! http://cassandra.nd.edu +! Prof. Edward Maginn +! Copyright (2013) University of Notre Dame du Lac +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see . +!******************************************************************************* + +SUBROUTINE Shell_GEMC_Particle_Transfer(box_in, box_out) + + !******************************************************************************** + ! + ! This subroutine performs particle swaps in a GEMC ensemble when core-shell units are present in the box + ! + ! + ! Called by + ! + ! gemc_driver + !********************************************************************************* + + USE Run_Variables + USE Random_Generators + USE Simulation_Properties + USE Energy_Routines + USE IO_Utilities + USE Fragment_Growth + USE Pair_Nrg_Routines + USE Volume + + IMPLICIT NONE + + INTEGER :: box_in, box_out, i, this_species, i_type + + INTEGER :: nmols_sorbate, nsorbate, is, ibox, im ,alive, which_anchor + INTEGER :: rand_igas + + INTEGER, ALLOCATABLE :: sorbate_id(:), frag_order(:) + + REAL(DP), ALLOCATABLE :: sorbate_x(:) + + REAL(DP) :: pick_species, delta_e_out, delta_e_out_pacc, e_bond_out, e_angle_out, e_dihed_out + REAL(DP) :: E_intra_vdw_out, E_intra_qq_out + REAL(DP) :: E_inter_vdw_out, E_inter_qq_out + REAL(DP) :: E_intra_vdw_igas, E_intra_qq_igas + REAL(DP) :: E_reciprocal_move_out, E_self_move_out, e_lrc_out, e_improper_out + REAL(DP) :: E_intra_vdw_in, E_intra_qq_in + REAL(DP) :: E_inter_vdw_in, E_inter_qq_in + REAL(DP) :: E_reciprocal_move_in, E_self_move_in, e_lrc_in, E_bond_in, E_improper_in + REAL(DP) :: E_angle_in, E_dihed_in, delta_e_in, delta_e_in_pacc, factor, P_forward, P_reverse, potw, CP_energy, P_seq + REAL(DP) :: lambda_for_build + LOGICAL :: inter_overlap, accept, accept_or_reject, cbmc_overlap, forward + LOGICAL :: intra_overlap + + REAL(DP) :: nrg_ring_frag_in, nrg_ring_frag_out + + TYPE(atom_class), ALLOCATABLE :: new_atom_list(:) + TYPE(molecule_class) :: new_molecule_list + TYPE(Energy_Class) :: energy_old_in, energy_old_out + + ! Variables added for l_pair_nrg and reciprocal k vector storage + + INTEGER :: position, box_number + + REAL(DP) :: time0, time1, pick_box + REAL(DP) :: E_vdw_move_alive, E_qq_move_alive + LOGICAL :: l_charge_in, l_charge_out, overlap + LOGICAL :: conv_Drude1, conv_Drude2 + INTEGER :: tot_mol_box_in, tot_mol_box_out + + potw = 1.0_DP + forward = .FALSE. + inter_overlap = .false. + cbmc_overlap = .false. + + P_seq = 1.0_DP + P_forward = 1.0_DP + P_reverse = 1.0_DP + + box_out = INT (rranf() * nbr_boxes) + 1 + + + pick_box = rranf() + + DO ibox = 1,nbr_boxes + IF (ibox == box_out) CYCLE + + IF (pick_box <= prob_swap_boxes(box_out,ibox)) EXIT + + END DO + + box_in = ibox + + + P_forward = P_forward * prob_swap_boxes(box_out,box_in) + P_reverse = P_reverse * prob_swap_boxes(box_in,box_out) + + ! increment total number of trials for each of the boxes + tot_trials(box_out) = tot_trials(box_out) + 1 + tot_trials(box_in) = tot_trials(box_in) + 1 + + ! Choose a species based on the overall mole fraction. + ! We will base this on the sorbate mole fraction + + IF (l_mol_frac_swap) THEN + + ! Choose a species based on the overall mole fraction. + ! We will base this on the sorbate mole fraction + ALLOCATE(sorbate_id(nspecies), sorbate_x(nspecies)) + + nmols_sorbate = 0 + nsorbate = 0 + sorbate_x(:) = 0.0_DP + + DO is = 1, nspecies + nsorbate = nsorbate + 1 + sorbate_id(nsorbate) = is + nmols_sorbate = nmols_sorbate + SUM(nmols(is,:)) + sorbate_x(nsorbate) = REAL(SUM(nmols(is,:)),DP) + END DO + + sorbate_x(:) = sorbate_x(:) / nmols_sorbate + + DO i = 2, nsorbate + sorbate_x(i) = sorbate_x(i) + sorbate_x(i-1) + END DO + + ! Now pick a species at random using Golden sampling + + pick_species = rranf() + + DO i = 1, nsorbate + IF (pick_species < sorbate_x(i)) EXIT + END DO + + this_species = sorbate_id(i) + + ELSE + + ! pick the species based on specified probability + + pick_species = rranf() + + DO is = 1, nspecies + IF (pick_species <= prob_swap_species(is)) EXIT + END DO + + this_species = is + + END IF + + + IF(nmols(this_species,box_out) == 0) THEN + ! IF(cpcollect) CALL Chempot(box_in) + RETURN + END IF + + ntrials(this_species,box_out)%deletion = ntrials(this_species,box_out)%deletion + 1 + ntrials(this_species,box_in)%insertion = ntrials(this_species,box_in)%insertion + 1 + ntrials(this_species,box_in)%cpcalc = ntrials(this_species,box_in)%cpcalc + 1 + + IF( box_list(box_out)%volume .lt. box_list(box_in)%volume ) forward = .TRUE. + + ! pick a molecule at random to delete + + im = INT(rranf() * nmols(this_species,box_out)) + 1 + + CALL Get_Index_Molecule(box_out,this_species,im,alive) + + CALL Save_Cartesian_Coordinates_Box(box_out,tot_mol_box_out) + CALL Save_Cartesian_Coordinates_Box(box_in,tot_mol_box_in) + + + l_charge_in = .FALSE. + l_charge_out = .FALSE. + IF((int_charge_sum_style(box_in) == charge_ewald .or. int_charge_sum_style(box_in) == charge_gaussian) .AND. has_charge(this_species)) l_charge_in = .TRUE. + IF((int_charge_sum_style(box_out) == charge_ewald .or. int_charge_sum_style(box_out) == charge_gaussian) .AND. has_charge(this_species)) l_charge_out = .TRUE. + + ! Let us first insert the molecule in other box and determine + ! if there is an overlap. Computationally this is more efficient + ! in the event that an overlap is detected, the move can be + ! immediately rejected. Note that this move chanegs the box identity + ! and coordinate related quantites of 'alive. When we calcualte + ! the energy of removing molecule from box_out we need to restore + ! + + molecule_list(alive,this_species)%which_box = box_in + + ! set deletion flag to false and call the CBMC growth routine + + cbmc_overlap = .FALSE. + + delta_e_in = 0.0_DP + + IF (species_list(this_species)%fragment .AND. species_list(this_species)%int_insert .NE. int_igas ) THEN + + del_Flag = .FALSE. + get_fragorder = .TRUE. + + ALLOCATE(frag_order(nfragments(this_species))) + ! frag_order = 0 + ! print *, 'frag_order',frag_order + lambda_for_build = molecule_list(alive,this_species)%cfc_lambda + !print *, species_list(this_species)%lcom + + IF (species_list(this_species)%lcom) THEN + CALL Build_Molecule(alive,this_species,box_in,frag_order,lambda_for_build, P_seq,P_forward, & + nrg_ring_frag_in,cbmc_overlap) + ELSE + CALL Build_Rigid_Fragment(alive,this_species,box_in,frag_order,lambda_for_build, P_seq,P_forward, & + nrg_ring_frag_in,cbmc_overlap) + END IF + ELSE + CALL New_Positions(box_in,alive,this_species,rand_igas) + END IF + + ! print *, frag_order + CALL Get_COM(alive,this_species) + CALL Compute_Max_COM_Distance(alive,this_species) + + CALL Fold_Molecule(alive,this_species,box_in) + + ! Note that the call to Build_Molecule returns cbmc_overlap flag to be + ! true not only in the case of an overlap but also when the weight + ! of a given fragment is zero. + + IF (.NOT. cbmc_overlap) THEN + + CALL compute_molecule_inter_overlap(alive,this_species,box_in, inter_overlap) + + END IF + + IF (inter_overlap .OR. cbmc_overlap) THEN + ! change the box identity + molecule_list(alive,this_species)%which_box = box_out + ! Revert the coordinates + + CALL Revert_Old_Cartesian_Coordinates(alive,this_species) + + ! It may so happen that only a few atoms have been assigned the exist flag to be true + ! during molecule building so make all the atom beads of alive as true. Also + ! reset the cfc_lambda + + atom_list(:,alive,this_species)%exist = .TRUE. + molecule_list(alive,this_species)%cfc_lambda = 1.0_DP + + RETURN + + END IF + + IF (int_vdw_sum_style(box_in) == vdw_cut_tail .or. int_vdw_sum_style(box_in) == born_cut_tail) THEN + nbeads_in(:) = nint_beads(:,box_in) + + DO i = 1, natoms(this_species) + i_type = nonbond_list(i,this_species)%atom_type_number + nint_beads(i_type,box_in) = nint_beads(i_type,box_in) + 1 + END DO + + END IF + + nmols(this_species,box_in) = nmols(this_species,box_in) + 1 + energy_old_in = energy(box_in) + + conv_drude1 = .TRUE. + + + IF (shell_mpm) CALL shell_relax(box_in, conv_drude1) + IF (conv_drude1) CALL Compute_Total_System_Energy(box_in, .TRUE., overlap) + + nmols(this_species,box_in) = nmols(this_species,box_in) - 1 + + ! Now compute the energy of the molecule assuming that it has been + ! deleted. Now if we are at this stage then we need to preserve + ! coordinates of the molecule inserted in other box. + + ALLOCATE(new_atom_list(natoms(this_species))) + + new_atom_list(:) = atom_list(:,alive,this_species) + new_molecule_list = molecule_list(alive,this_species) + + + ! Now we will swap the particle in another box and compute the energy + ! For this, we will change the box identity for this molecule + ! Now revert to the original coordinates + + CALL Revert_Old_Cartesian_Coordinates(alive,this_species) + molecule_list(alive,this_species)%which_box = box_out + + + delta_e_out = 0.0_DP + + + IF ( species_list(this_species)%fragment .AND. species_list(this_species)%int_insert .NE. int_igas) THEN + ! Note that we need to use the same fragment order in which we inserted the molecule + ! above. so in this case frag_order becomes input to the routine. We obtain + ! P_delete via this call. Note that, cbmc_overlap should be false as we are + ! dealing with an existing molecule. + del_Flag = .TRUE. + get_fragorder = .FALSE. + + IF (species_list(this_species)%lcom) THEN + CALL Build_Molecule(alive,this_species,box_out,frag_order,lambda_for_build, P_seq,P_reverse, & + nrg_ring_frag_out,cbmc_overlap) + ELSE + CALL Build_Rigid_Fragment(alive,this_species,box_out,frag_order,lambda_for_build, & + P_seq, P_reverse, nrg_ring_frag_out,cbmc_overlap) + END IF + +! DEALLOCATE(frag_order) + + IF (cbmc_overlap) THEN + ! energy overlap was found in the old configuration + + atom_list(1:natoms(this_species),alive,this_species)%exist = .TRUE. + CALL Revert_Old_Cartesian_Coordinates(alive,this_species) + + molecule_list(alive,this_species)%cfc_lambda = 1.0_DP + + END IF + + + END IF + + molecule_list(alive,this_species)%which_box = box_in + + IF (int_vdw_sum_style(box_out) == vdw_cut_tail .or. int_vdw_sum_style(box_out) == born_cut_tail) THEN + + nbeads_out(:) = nint_beads(:,box_out) + + DO i = 1, natoms(this_species) + i_type = nonbond_list(i,this_species)%atom_type_number + nint_beads(i_type,box_out) = nint_beads(i_type,box_out) - 1 + END DO + + END IF + + nmols(this_species,box_out) = nmols(this_species,box_out) - 1 + energy_old_out = energy(box_out) + conv_drude2 = .TRUE. + + IF(shell_mpm) CALL shell_relax(box_out, conv_drude2) + IF (conv_drude2) CALL Compute_Total_System_Energy(box_out, .TRUE., overlap) + + nmols(this_species,box_out) = nmols(this_species,box_out) + 1 + + ! Define a factor that will be used to accept or reject the move. Note that + ! the change in energy of box_out is actually negative of delta_e_out calculated + ! above + + !factor = beta(box_in)*delta_e_in_pacc - beta(box_out)*delta_e_out_pacc + + !factor = beta(box_in)*(energy(box_in)%total - energy_old_in%total - e_angle_in - nrg_ring_frag_in) & + ! & + beta(box_out)*(energy(box_out)%total - energy_old_out%total - e_angle_out - nrg_ring_frag_out) + + + factor = beta(box_in)*(energy(box_in)%total - energy_old_in%total ) + beta(box_out)*(energy(box_out)%total - energy_old_out%total) + + !IF(dabs(e_angle_in)> 1d-5 .or. dabs(nrg_ring_frag_in)>1d-5 .or. dabs(e_angle_out) > 1d-5 .or. dabs(nrg_ring_frag_out)>1d-5) THEN + ! print *, 'nrg_ring_frag_in .NE. 0.0', nrg_ring_frag_in, e_angle_in,e_angle_out, nrg_ring_frag_out + !END IF + + factor = factor - DLOG((box_list(box_in)%volume * REAL(nmols(this_species,box_out),DP))/ & + (box_list(box_out)%volume * REAL(nmols(this_species,box_in) + 1, DP))) + + factor = factor + DLOG(P_forward) - DLOG(P_reverse) + + accept = accept_or_reject(factor) + + IF (accept .and. conv_drude1 .and. conv_Drude2) THEN + + nmols(this_species,box_in) = nmols(this_species,box_in) + 1 + nmols(this_species,box_out) = nmols(this_species,box_out) - 1 + + nsuccess(this_species,box_in)%insertion = nsuccess(this_species,box_in)%insertion + 1 + nsuccess(this_species,box_out)%deletion = nsuccess(this_species,box_out)%deletion + 1 + + DO i = 1,natoms(this_species) + atom_list(i,alive,this_species) = new_atom_list(i) + END DO + molecule_list(alive,this_species) = new_molecule_list + + CALL Fold_Molecule(alive,this_species,box_in) + + ! Restore the coordinates of the molecule due to successful insertion + + CALL Get_Internal_Coordinates(alive,this_species) + + + ELSE + + molecule_list(alive,this_species)%which_box = box_out + IF ( int_vdw_sum_style(box_in) == vdw_cut_tail .or. int_vdw_sum_style(box_in) == born_cut_tail ) THEN + + nint_beads(:,box_in) = nbeads_in(:) + + END IF + + IF ( int_vdw_sum_style(box_out) == vdw_cut_tail .or. int_vdw_sum_style(box_in) == born_cut_tail ) THEN + + nint_beads(:,box_out) = nbeads_out(:) + + END IF + + CALL Reset_Cartesian_Coordinates_Box(box_in) + CALL Reset_Cartesian_Coordinates_Box(box_out) + + energy(box_out) = energy_old_out + energy(box_in) = energy_old_in + + END IF + + + IF (ALLOCATED(sorbate_id)) DEALLOCATE(sorbate_id) + IF (ALLOCATED(sorbate_x)) DEALLOCATE(sorbate_x) + DEALLOCATE(new_atom_list) + +END SUBROUTINE Shell_GEMC_Particle_Transfer + + diff --git a/Src/shell_relax.f90 b/Src/shell_relax.f90 new file mode 100644 index 00000000..83513b56 --- /dev/null +++ b/Src/shell_relax.f90 @@ -0,0 +1,101 @@ +SUBROUTINE shell_relax (this_box, conv_drude) + +USE Type_Definitions +USE Run_Variables +USE Random_Generators +USE Simulation_Properties +USE Energy_Routines +!USE Pair_Nrg_Routines + +IMPLICIT NONE + +INTEGER, INTENT(IN) :: this_box +LOGICAL, INTENT(OUT) :: conv_drude + +INTEGER :: is, im, ia, alive +REAL(DP) :: deltaL, xdold, ydold, zdold +INTEGER :: iter_drude, this_locate +REAL(DP) :: lambda_it + + +lambda_it = 1.0 +iter_drude = 0 +conv_drude = .FALSE. + +DO WHILE ( (.NOT. conv_drude) .and. iter_drude .LE. 100) !Perform 100 iterations + +iter_drude = iter_drude + 1 + +CALL Compute_drude_cos_sin_sum(this_box) + +DO is = 1, nspecies + DO ia = 1, natoms(is) + IF(atom_list(ia,1,is)%drude_type) THEN + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP PRIVATE(im, this_locate) & + !$OMP SCHEDULE(STATIC) + DO im = 1, nmolecules(is) + this_locate = locate(im,is) + IF( .NOT. molecule_list(this_locate,is)%live) CYCLE + IF( molecule_list(this_locate,is)%which_box /= this_box ) CYCLE + CALL Compute_molecule_drude_force(is,this_locate,ia,this_box, drudeFx(is,this_locate,ia), drudeFy(is,this_locate,ia), drudeFz(is,this_locate,ia) ) + END DO + !$OMP END PARALLEL DO + END IF + END DO +END DO + +!if(iter_drude > 20) lambda_it = 0.1 !if the EM cannot converge in 20 steps then use a smaller step size. This is usually due to close intermolecular distance or soft harmonic spring. + + +conv_drude = .TRUE. + +DO is = 1, nspecies + DO ia = 1, natoms(is) + IF(atom_list(ia,1,is)%drude_type) THEN + !$OMP PARALLEL DO DEFAULT(SHARED) & + !$OMP PRIVATE(im, this_locate,deltaL) & + !$OMP PRIVATE(xdold, ydold, zdold) & + !$OMP SCHEDULE(STATIC) & + !$OMP REDUCTION(.AND.:conv_drude) + DO im = 1, nmolecules(is) !nmols(is,this_box) !nmolecules(is) + this_locate = locate(im,is) + IF( .NOT. molecule_list(this_locate,is)%live) CYCLE + IF( molecule_list(this_locate,is)%which_box /= this_box ) CYCLE + + deltaL = 0.0 + xdold = atom_list(ia,this_locate,is)%rxp + ydold = atom_list(ia,this_locate,is)%ryp + zdold = atom_list(ia,this_locate,is)%rzp + + atom_list(ia,this_locate,is)%rxp = xdold + drudeFx(is,this_locate,ia)/nonbond_list(ia,is)%pol_alpha*lambda_it + atom_list(ia,this_locate,is)%ryp = ydold + drudeFy(is,this_locate,ia)/nonbond_list(ia,is)%pol_alpha*lambda_it + atom_list(ia,this_locate,is)%rzp = zdold + drudeFz(is,this_locate,ia)/nonbond_list(ia,is)%pol_alpha*lambda_it + + deltaL = (atom_list(ia,this_locate,is)%rxp - xdold)*(atom_list(ia,this_locate,is)%rxp - xdold) +(atom_list(ia,this_locate,is)%ryp - ydold)*(atom_list(ia,this_locate,is)%ryp - ydold) + (atom_list(ia,this_locate,is)%rzp - zdold)*(atom_list(ia,this_locate,is)%rzp - zdold) + + if(sqrt(deltaL) > 1d-3 .or. maxval(drudeFx) > 50 .or. maxval(drudeFy) > 50.0 .or. maxval(drudeFz) > 50.0 ) then !convergence check + conv_drude = .FALSE. + end if + END DO + !$OMP END PARALLEL DO + END IF + END DO +END DO + +END DO + +IF(.NOT. conv_drude) THEN +check_err = check_err + 1 !accumulate unconvergence counter +print *, 'EM did not converge in 100 iterations' +END IF + +END SUBROUTINE shell_relax + + + + + + + + diff --git a/Src/type_definitions.f90 b/Src/type_definitions.f90 index aa6e4be4..17b198fa 100644 --- a/Src/type_definitions.f90 +++ b/Src/type_definitions.f90 @@ -200,6 +200,7 @@ MODULE Type_Definitions REAL(DP) :: rxp_nls, ryp_nls, rzp_nls ! The starting positions for the neighbor list REAL(DP) :: rxp_old, ryp_old, rzp_old LOGICAL :: exist + LOGICAL :: drude_type END TYPE Atom_Class !**************************************************************************** @@ -229,6 +230,10 @@ MODULE Type_Definitions LOGICAL :: ring_atom + REAL(DP) :: g_alpha + REAL(DP) :: pol_alpha + + END TYPE Nonbond_Class !**************************************************************************** diff --git a/Src/volume_change.f90 b/Src/volume_change.f90 index 6f865d78..21ae71fb 100644 --- a/Src/volume_change.f90 +++ b/Src/volume_change.f90 @@ -84,7 +84,7 @@ SUBROUTINE Volume_Change(this_box,this_step) REAL(DP) :: pres_id REAL(DP) :: W_vol_vdw, W_vol_qq - LOGICAL :: accept, allocation_cos_sin, overlap, xz_change, accept_or_reject + LOGICAL :: accept, allocation_cos_sin, overlap, xz_change, accept_or_reject, conv_drude TYPE(Box_Class) :: box_list_old @@ -178,7 +178,7 @@ SUBROUTINE Volume_Change(this_box,this_step) END IF - IF (int_charge_sum_style(this_box) == charge_ewald) THEN + IF (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian ) THEN ! store the cos_mol and sin_mol arrays @@ -244,7 +244,7 @@ SUBROUTINE Volume_Change(this_box,this_step) rcut_max_old = rcut_max(this_box) - IF( int_charge_sum_style(this_box) == charge_ewald ) THEN + IF( int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian ) THEN ! alpha_ewald_old = alpha_ewald(this_box) h_ewald_cut_old = h_ewald_cut(this_box) @@ -263,7 +263,7 @@ SUBROUTINE Volume_Change(this_box,this_step) rcut_max(this_box) = rcut_vdw(this_box) - IF ( int_charge_sum_style(this_box) == charge_ewald) THEN + IF ( int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian) THEN ! alpha_ewald(this_box) = ewald_p_sqrt(this_box) / rcut_coul(this_box) h_ewald_cut(this_box) = 2.0_DP * ewald_p(this_box) / rcut_coul(this_box) @@ -361,7 +361,7 @@ SUBROUTINE Volume_Change(this_box,this_step) energy_old = energy(this_box) - IF (int_charge_sum_style(this_box) == charge_ewald) THEN + IF (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian) THEN ! The calculation of Ewald reciprocal energy is tricky. The k vectors for this_box will change along with ! total number of k vectors for the box. It may so happen that the resultant k vectors may exceed the @@ -450,14 +450,16 @@ SUBROUTINE Volume_Change(this_box,this_step) END IF + + conv_drude = .TRUE. + IF(shell_mpm) THEN + CALL shell_relax(this_box, conv_drude) + IF (conv_drude) CALL Compute_Total_System_Energy(this_box,.TRUE.,overlap) + ELSE + CALL Compute_Total_System_Energy(this_box,.TRUE.,overlap) + END IF - CALL Compute_Total_System_Energy(this_box,.TRUE.,overlap) - - - - - - + IF (overlap) THEN ! reject move CALL Reset_Coords @@ -468,8 +470,7 @@ SUBROUTINE Volume_Change(this_box,this_step) delta_e = energy(this_box)%total - energy_old%total - ! based on the energy, calculate the acceptance ratio - + ! based on the energy, calculate the acceptance ratio ln_pacc = beta(this_box) * delta_e & + beta(this_box) * pressure(this_box) * delta_volume & @@ -477,15 +478,14 @@ SUBROUTINE Volume_Change(this_box,this_step) accept = accept_or_reject(ln_pacc) - IF ( accept ) THEN - + IF ( accept .and. conv_drude) THEN nvol_success(this_box) = nvol_success(this_box) + 1 ivol_success(this_box) = ivol_success(this_box) + 1 ! energy, positions and box dimensions are already updated - IF (int_charge_sum_style(this_box) == charge_ewald) THEN + IF (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian) THEN ! put the sin_sum and cos_sum for other boxes ! Note that hx,hy,hz and Cn are updated when the total energy routine is called above. So ! there is no need to update these arrays. @@ -655,7 +655,7 @@ SUBROUTINE Reset_Coords rcut_max_old = rcut_max(this_box) - IF( int_charge_sum_style(this_box) == charge_ewald ) THEN + IF( int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian ) THEN ! alpha_ewald(this_box) = alpha_ewald_old h_ewald_cut(this_box) = h_ewald_cut_old @@ -665,7 +665,7 @@ SUBROUTINE Reset_Coords END IF - IF (int_charge_sum_style(this_box) == charge_ewald) THEN + IF (int_charge_sum_style(this_box) == charge_ewald .or. int_charge_sum_style(this_box) == charge_gaussian ) THEN ! reset the terms related to Ewald reciprocal energy ! reset the total number of kvectors diff --git a/Src/write_properties.f90 b/Src/write_properties.f90 index eaf14554..ab5a7a87 100644 --- a/Src/write_properties.f90 +++ b/Src/write_properties.f90 @@ -172,7 +172,7 @@ SUBROUTINE Write_Properties_Buffer(file_number) USE Simulation_Properties - INTEGER :: file_number, ii, is, is_dens, is_cp, is_lambda, total_frac + INTEGER :: file_number, ii, is, is_dens, is_cp, is_lambda, total_frac, ibox, jbox REAL(DP),DIMENSION(:), ALLOCATABLE :: write_buff CHARACTER(FILENAME_LEN) :: prop_written @@ -195,6 +195,25 @@ SUBROUTINE Write_Properties_Buffer(file_number) is_cp = 1 is_lambda = 1 + IF (int_sim_type == sim_gemc .and. this_box == 1) THEN + + constant_vol = .TRUE. + + p_plus_1 = .TRUE. + p_calc = .FALSE. + DO WHILE(.NOT. p_calc) + CALL shell_gemc_nvt_volume(ibox,jbox) + END DO + p_plus_1 = .FALSE. + p_calc = .FALSE. + DO WHILE(.NOT. p_calc) + CALL shell_gemc_nvt_volume(ibox,jbox) + END DO + + constant_vol = .FALSE. + + END IF + DO WHILE ( ii <= prop_per_file(file_number,this_box)) prop_written = prop_output(ii,file_number,this_box) @@ -245,7 +264,7 @@ SUBROUTINE Write_Properties_Buffer(file_number) P_inst(this_box) = ((Pressure_tensor(1,1,this_box) + Pressure_tensor(2,2,this_box) + & Pressure_tensor(3,3,this_box)) / 3.0_DP) * atomic_to_bar - IF(int_vdw_sum_style(this_box) == vdw_cut_tail) THEN + IF(int_vdw_sum_style(this_box) == vdw_cut_tail .or. int_vdw_sum_style(this_box) == born_cut_tail) THEN P_inst(this_box) = P_inst(this_box) + ((virial(this_box)%lrc / box_list(this_box)%volume) * atomic_to_bar) END IF @@ -253,6 +272,14 @@ SUBROUTINE Write_Properties_Buffer(file_number) write_buff(ii+1) = P_ideal(this_box) + P_inst(this_box) + ELSE IF (prop_written == 'Pressure1') THEN + + write_buff(ii+1) = P_inst_plus(this_box) + + ELSE IF (prop_written == 'Pressure2') THEN + + write_buff(ii+1) = P_inst_minus(this_box) + ELSE IF (prop_written == 'Volume') THEN IF (block_average) THEN