From 3ed860207c25034330aabb5e04d8020f95c29de8 Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Tue, 12 Oct 2021 15:27:24 -0500 Subject: [PATCH 01/16] Added comments to explain eventual electorstatic fix --- src/Ewald.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/Ewald.cpp b/src/Ewald.cpp index 80a879249..b7202bf43 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -909,12 +909,18 @@ double Ewald::MolCorrection(uint molIndex, uint box) const continue; } for (uint j = i + 1; j < atomSize; j++) { - currentAxes.InRcut(distSq, virComponents, currentCoords, + if(currentAxes.InRcut(distSq, virComponents, currentCoords, start + i, start + j, box); dist = sqrt(distSq); correction += (thisKind.AtomCharge(i) * thisKind.AtomCharge(j) * erf(ff.alpha[box] * dist) / dist); } + // within cutoff + // CalcCoul = qiqj/r(erfc(alpha*r)+erf(alpha*r)) + + // outside cutoff + // erf(alpha*r)*qiqj/r + // CalcCoul = erf } return -1.0 * num::qqFact * correction * lambdaCoef * lambdaCoef; From 3fdbdbfda60b6e624ede16fe2f0ef94a2127de4a Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Mon, 22 Nov 2021 13:52:26 -0600 Subject: [PATCH 02/16] wrote scripts to automate integration testing --- test/Run_Examples.py | 109 +++++++++++++++++++++++++++++++++++++++++ test/Setup_Examples.sh | 44 +++++++++++++++++ 2 files changed, 153 insertions(+) create mode 100644 test/Run_Examples.py create mode 100755 test/Setup_Examples.sh diff --git a/test/Run_Examples.py b/test/Run_Examples.py new file mode 100644 index 000000000..8601a4070 --- /dev/null +++ b/test/Run_Examples.py @@ -0,0 +1,109 @@ +#!/usr/bin/python + +def substring_after(s, delim): + return s.partition(delim)[2] + +import os +import glob +import sys +import subprocess +import pandas as pd + +CPU_binaries_dict = {} +GPU_binaries_dict = {} + +os.chdir("new_feature_binaries") +pathToDir = os.getcwd() +CPU_binaries_new_feature = sorted(glob.glob('*_CPU_*'), key=os.path.getmtime) +GPU_binaries_new_feature = sorted(glob.glob('*_GPU_*'), key=os.path.getmtime) + +for binary in CPU_binaries_new_feature: + CPU_binaries_dict[substring_after(binary, "GOMC_CPU_")+"_new"] = os.path.join(pathToDir, binary) + +for binary in GPU_binaries_new_feature: + GPU_binaries_dict[substring_after(binary, "GOMC_GPU_")+"_new"] = os.path.join(pathToDir, binary) + +os.chdir("../ref_binaries") +pathToDir = os.getcwd() +CPU_binaries_ref = sorted(glob.glob('*_CPU_*'), key=os.path.getmtime) +GPU_binaries_ref = sorted(glob.glob('*_GPU_*'), key=os.path.getmtime) + +for binary in CPU_binaries_ref: + CPU_binaries_dict[substring_after(binary, "GOMC_CPU_")+"_ref"] = os.path.join(pathToDir, binary) + +for binary in GPU_binaries_ref: + GPU_binaries_dict[substring_after(binary, "GOMC_GPU_")+"_ref"] = os.path.join(pathToDir, binary) + + +print("CPU Binaries Dict", CPU_binaries_dict) +print("GPU Binaries Dict", GPU_binaries_dict) + +os.chdir("../integration") + +confFileName = "in.conf" + +Log_Template_file = open("IntegrationTest.log", 'w') +print("opened Log") + +listOfTests = [] + +# traverse root directory, and list directories as dirs and files as files +# traverse root directory, and list directories as dirs and files as files +for root, dirs, files in os.walk("."): + path = root.split(os.sep) +# print((len(path) - 1) * '---', os.path.basename(root)) + for file in files: +# print(len(path) * '---', file) + if file==confFileName: + newOrRef = "" + if "new_feature_cpu" in path: + newOrRef = "_new" + elif "ref_cpu" in path: + newOrRef = "_ref" + + run = False + if "NVT"+newOrRef in CPU_binaries_dict and 'NVT' in path and 'GEMC_NVT' not in path: + print("Call GOMC") + command = CPU_binaries_dict["NVT"+newOrRef],os.path.abspath(root),os.path.basename(root) + print(CPU_binaries_dict["NVT"+newOrRef],os.path.abspath(root),os.path.basename(root)) + listOfTests.append(command) + run = True + elif "NPT"+newOrRef in CPU_binaries_dict and 'NPT' in path and 'GEMC_NPT' not in path: + print("Call GOMC") + print(CPU_binaries_dict["NPT"+newOrRef],os.path.abspath(root),os.path.basename(root)) + command = CPU_binaries_dict["NPT"+newOrRef],os.path.abspath(root),os.path.basename(root) + listOfTests.append(command) + run = True + elif "GCMC"+newOrRef in CPU_binaries_dict and 'GCMC' in path: + print("Call GOMC") + print(CPU_binaries_dict["GCMC"+newOrRef],os.path.abspath(root),os.path.basename(root)) + command = CPU_binaries_dict["GCMC"+newOrRef],os.path.abspath(root),os.path.basename(root) + listOfTests.append(command) + run = True + elif "GEMC"+newOrRef in CPU_binaries_dict and 'GEMC_NVT' in path: + print("Call GOMC") + print(CPU_binaries_dict["GEMC"+newOrRef],os.path.abspath(root),os.path.basename(root)) + command = CPU_binaries_dict["GEMC"+newOrRef],os.path.abspath(root),os.path.basename(root) + listOfTests.append(command) + run = True + elif "GEMC"+newOrRef in CPU_binaries_dict and 'GEMC_NPT' in path: + print("Call GOMC") + print(CPU_binaries_dict["GEMC"+newOrRef],os.path.abspath(root),os.path.basename(root)) + command = CPU_binaries_dict["GEMC"+newOrRef],os.path.abspath(root),os.path.basename(root) + listOfTests.append(command) + run = True + else: + run = False + +print(listOfTests) +# Create the pandas DataFrame +df = pd.DataFrame(listOfTests, columns = ['PathToBinary', 'PathToExample', 'Example']) +df = df.sort_values(by=['Example']) +print(df) +# if(run): +# exec_GOMC_run_command = subprocess.Popen(command, shell=True, stderr=subprocess.STDOUT) +# write_log_data = f'Waiting for GOMC Example {command} to finish.' +# Log_Template_file.write(str(write_log_data)) +# print(str(write_log_data)) +# GOMC_pid_status = os.waitpid(exec_GOMC_run_command.pid, os.WSTOPPED) # pauses python until box 0 sim done +# write_log_data = f'The GOMC Example {command} is finished.' diff --git a/test/Setup_Examples.sh b/test/Setup_Examples.sh new file mode 100755 index 000000000..e3182e863 --- /dev/null +++ b/test/Setup_Examples.sh @@ -0,0 +1,44 @@ +#!/bin/bash +git clone https://github.com/GOMC-WSU/GOMC_Examples.git +mkdir integration +cd integration + +mkdir new_feature_cpu +cp -frd GOMC_Examples new_feature_cpu + +mkdir ref_cpu +cp -frd GOMC_Examples ref_cpu + +mkdir new_feature_gpu +cp -frd GOMC_Examples new_feature_gpu + +mkdir ref_gpu +cp -frd GOMC_Examples ref_gpu + + +cd ../.. +branch=$(git symbolic-ref HEAD | sed -e 's,.*/\(.*\),\1,') +echo "Building $branch binaries" + +./metamake.sh -g NVT +mkdir -p test/new_feature_binaries +cp -frdp ./bin/* test/new_feature_binaries + + +echo "$branch" +if [ $branch == "development" ]; then + echo "I am on development; checking out main" + git checkout main +else + echo "I am on $branch; checking out development" + git checkout development +fi + +rm -frd bin +branch=$(git symbolic-ref HEAD | sed -e 's,.*/\(.*\),\1,') +echo "Building $branch binaries" +./metamake.sh -g NVT +mkdir -p test/ref_binaries +cp -frd ./bin/* test/ref_binaries +cd test + From 6557396e04166d3d41a88a0b806d1035add2a33d Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Mon, 22 Nov 2021 18:24:18 -0600 Subject: [PATCH 03/16] Set --- test/Setup_Examples.sh | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/test/Setup_Examples.sh b/test/Setup_Examples.sh index e3182e863..d5859e467 100755 --- a/test/Setup_Examples.sh +++ b/test/Setup_Examples.sh @@ -3,40 +3,40 @@ git clone https://github.com/GOMC-WSU/GOMC_Examples.git mkdir integration cd integration -mkdir new_feature_cpu -cp -frd GOMC_Examples new_feature_cpu +mkdir new_cpu +cp -frd GOMC_Examples new_cpu mkdir ref_cpu cp -frd GOMC_Examples ref_cpu -mkdir new_feature_gpu -cp -frd GOMC_Examples new_feature_gpu +mkdir new_gpu +cp -frd GOMC_Examples new_gpu mkdir ref_gpu cp -frd GOMC_Examples ref_gpu cd ../.. -branch=$(git symbolic-ref HEAD | sed -e 's,.*/\(.*\),\1,') -echo "Building $branch binaries" +startingBranch=$(git symbolic-ref HEAD | sed -e 's,.*/\(.*\),\1,') +echo "Building $startingBranch binaries" ./metamake.sh -g NVT -mkdir -p test/new_feature_binaries -cp -frdp ./bin/* test/new_feature_binaries +mkdir -p test/new_binaries +cp -frdp ./bin/* test/new_binaries -echo "$branch" -if [ $branch == "development" ]; then +echo "$startingBranch" +if [ $startingBranch == "development" ]; then echo "I am on development; checking out main" git checkout main else - echo "I am on $branch; checking out development" + echo "I am on $startingBranch; checking out development" git checkout development fi rm -frd bin -branch=$(git symbolic-ref HEAD | sed -e 's,.*/\(.*\),\1,') -echo "Building $branch binaries" +refBranch=$(git symbolic-ref HEAD | sed -e 's,.*/\(.*\),\1,') +echo "Building $refBbranch binaries" ./metamake.sh -g NVT mkdir -p test/ref_binaries cp -frd ./bin/* test/ref_binaries From e09bf883902b30572224a89d78845f91443caa92 Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Mon, 22 Nov 2021 18:27:03 -0600 Subject: [PATCH 04/16] update Run --- test/Run_Examples.py | 104 ++++++++++++++++++++++++------------------- 1 file changed, 59 insertions(+), 45 deletions(-) diff --git a/test/Run_Examples.py b/test/Run_Examples.py index 8601a4070..fc4aa1b05 100644 --- a/test/Run_Examples.py +++ b/test/Run_Examples.py @@ -6,37 +6,33 @@ def substring_after(s, delim): import os import glob import sys -import subprocess import pandas as pd - -CPU_binaries_dict = {} +import time +import datetime +import subprocess +from subprocess import Popen, PIPE, STDOUT +binaries_dict = {} GPU_binaries_dict = {} os.chdir("new_feature_binaries") pathToDir = os.getcwd() -CPU_binaries_new_feature = sorted(glob.glob('*_CPU_*'), key=os.path.getmtime) -GPU_binaries_new_feature = sorted(glob.glob('*_GPU_*'), key=os.path.getmtime) - -for binary in CPU_binaries_new_feature: - CPU_binaries_dict[substring_after(binary, "GOMC_CPU_")+"_new"] = os.path.join(pathToDir, binary) - -for binary in GPU_binaries_new_feature: - GPU_binaries_dict[substring_after(binary, "GOMC_GPU_")+"_new"] = os.path.join(pathToDir, binary) +binaries_cpu_new = sorted(glob.glob('GOMC_CPU_*'), key=os.path.getmtime) +binaries_gpu_new = sorted(glob.glob('GOMC_GPU_*'), key=os.path.getmtime) +for binary in binaries_cpu_new: + binaries_dict[substring_after(binary, "GOMC_")+"_new"] = os.path.join(pathToDir, binary) +for binary in binaries_gpu_new: + binaries_dict[substring_after(binary, "GOMC_")+"_new"] = os.path.join(pathToDir, binary) os.chdir("../ref_binaries") pathToDir = os.getcwd() -CPU_binaries_ref = sorted(glob.glob('*_CPU_*'), key=os.path.getmtime) -GPU_binaries_ref = sorted(glob.glob('*_GPU_*'), key=os.path.getmtime) - -for binary in CPU_binaries_ref: - CPU_binaries_dict[substring_after(binary, "GOMC_CPU_")+"_ref"] = os.path.join(pathToDir, binary) - -for binary in GPU_binaries_ref: - GPU_binaries_dict[substring_after(binary, "GOMC_GPU_")+"_ref"] = os.path.join(pathToDir, binary) - +binaries_cpu_ref = sorted(glob.glob('GOMC_CPU_*'), key=os.path.getmtime) +binaries_gpu_ref = sorted(glob.glob('GOMC_GPU_*'), key=os.path.getmtime) +for binary in binaries_cpu_ref: + binaries_dict[substring_after(binary, "GOMC_")+"_ref"] = os.path.join(pathToDir, binary) +for binary in binaries_gpu_ref: + binaries_dict[substring_after(binary, "GOMC_")+"_ref"] = os.path.join(pathToDir, binary) -print("CPU Binaries Dict", CPU_binaries_dict) -print("GPU Binaries Dict", GPU_binaries_dict) +print("Binaries Dict", binaries_dict) os.chdir("../integration") @@ -62,34 +58,34 @@ def substring_after(s, delim): newOrRef = "_ref" run = False - if "NVT"+newOrRef in CPU_binaries_dict and 'NVT' in path and 'GEMC_NVT' not in path: + if cpuOrGpu+"NVT"+newOrRef in binaries_dict and 'NVT' in path and 'GEMC_NVT' not in path: print("Call GOMC") - command = CPU_binaries_dict["NVT"+newOrRef],os.path.abspath(root),os.path.basename(root) - print(CPU_binaries_dict["NVT"+newOrRef],os.path.abspath(root),os.path.basename(root)) + command = binaries_dict[cpuOrGpu+"NVT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NVT"+newOrRef,os.path.basename(root) + print(binaries_dict[cpuOrGpu+"NVT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NVT"+newOrRef,os.path.basename(root)) listOfTests.append(command) run = True - elif "NPT"+newOrRef in CPU_binaries_dict and 'NPT' in path and 'GEMC_NPT' not in path: + elif cpuOrGpu+"NPT"+newOrRef in binaries_dict and 'NPT' in path and 'GEMC_NPT' not in path: print("Call GOMC") - print(CPU_binaries_dict["NPT"+newOrRef],os.path.abspath(root),os.path.basename(root)) - command = CPU_binaries_dict["NPT"+newOrRef],os.path.abspath(root),os.path.basename(root) + print(binaries_dict[cpuOrGpu+"NPT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NPT"+newOrRef,os.path.basename(root)) + command = binaries_dict[cpuOrGpu+"NPT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NPT"+newOrRef,os.path.basename(root) listOfTests.append(command) run = True - elif "GCMC"+newOrRef in CPU_binaries_dict and 'GCMC' in path: + elif cpuOrGpu+"GCMC"+newOrRef in binaries_dict and 'GCMC' in path: print("Call GOMC") - print(CPU_binaries_dict["GCMC"+newOrRef],os.path.abspath(root),os.path.basename(root)) - command = CPU_binaries_dict["GCMC"+newOrRef],os.path.abspath(root),os.path.basename(root) + print(binaries_dict[cpuOrGpu+"GCMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GCMC"+newOrRef,os.path.basename(root)) + command = binaries_dict[cpuOrGpu+"GCMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GCMC"+newOrRef,os.path.basename(root) listOfTests.append(command) run = True - elif "GEMC"+newOrRef in CPU_binaries_dict and 'GEMC_NVT' in path: + elif cpuOrGpu+"GEMC"+newOrRef in binaries_dict and 'GEMC_NVT' in path: print("Call GOMC") - print(CPU_binaries_dict["GEMC"+newOrRef],os.path.abspath(root),os.path.basename(root)) - command = CPU_binaries_dict["GEMC"+newOrRef],os.path.abspath(root),os.path.basename(root) + print(binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GEMC_NVT"+newOrRef,os.path.basename(root)) + command = binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GEMC_NVT"+newOrRef,os.path.basename(root) listOfTests.append(command) run = True - elif "GEMC"+newOrRef in CPU_binaries_dict and 'GEMC_NPT' in path: + elif cpuOrGpu+"GEMC"+newOrRef in binaries_dict and 'GEMC_NPT' in path: print("Call GOMC") - print(CPU_binaries_dict["GEMC"+newOrRef],os.path.abspath(root),os.path.basename(root)) - command = CPU_binaries_dict["GEMC"+newOrRef],os.path.abspath(root),os.path.basename(root) + print(binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GEMC_NPT"+newOrRef,os.path.basename(root)) + command = binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GEMC_NPT"+newOrRef,os.path.basename(root) listOfTests.append(command) run = True else: @@ -97,13 +93,31 @@ def substring_after(s, delim): print(listOfTests) # Create the pandas DataFrame -df = pd.DataFrame(listOfTests, columns = ['PathToBinary', 'PathToExample', 'Example']) +df = pd.DataFrame(listOfTests, columns = ['PathToBinary', 'PathToExample', 'Binary', 'Example']) df = df.sort_values(by=['Example']) print(df) -# if(run): -# exec_GOMC_run_command = subprocess.Popen(command, shell=True, stderr=subprocess.STDOUT) -# write_log_data = f'Waiting for GOMC Example {command} to finish.' -# Log_Template_file.write(str(write_log_data)) -# print(str(write_log_data)) -# GOMC_pid_status = os.waitpid(exec_GOMC_run_command.pid, os.WSTOPPED) # pauses python until box 0 sim done -# write_log_data = f'The GOMC Example {command} is finished.' +for index, row in df.iterrows(): + print(index) + print(row) +""" + os.chdir(row['PathToExample']) + write_log_data = "Changing directory to {}\n".format(row['PathToExample']) + Log_Template_file.write(str(write_log_data)) + command = (row['PathToBinary'] + " in.conf > out.log") + write_log_data = "Issuing command: {}\n".format(command) + Log_Template_file.write(str(write_log_data)) + start = time.time() + exec_GOMC_run_command = subprocess.Popen(command, shell=True, stderr=STDOUT) + write_log_data = "Waiting for GOMC Example {} {} to finish.\n".format(row['Binary'],row['Example']) + print(str(write_log_data)) + Log_Template_file.write(str(write_log_data)) + GOMC_pid_status = os.waitpid(exec_GOMC_run_command.pid, os.WSTOPPED) # pauses python until box 0 sim done + end = time.time() + write_log_data = "Elapsed time: {}.\n".format(datetime.timedelta(seconds=end-start)) + Log_Template_file.write(str(write_log_data)) + print(str(write_log_data)) + write_log_data = "The GOMC Example {} {} has finished.\n".format(row['Binary'],row['Example']) + print(str(write_log_data)) + Log_Template_file.write(str(write_log_data)) + Log_Template_file.flush() +""" From 047dc3cf40d6b8684ee07d1933707de32ccea6e2 Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Mon, 22 Nov 2021 18:57:36 -0600 Subject: [PATCH 05/16] change ewald back --- src/Ewald.cpp | 390 +++++++++++++++----------------------------------- 1 file changed, 114 insertions(+), 276 deletions(-) diff --git a/src/Ewald.cpp b/src/Ewald.cpp index 8d19fa151..6895148ed 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -24,7 +24,6 @@ along with this program, also can be found at . #include "CalculateForceCUDAKernel.cuh" #include "ConstantDefinitionsCUDAKernel.cuh" #endif -#include "GOMCEventsProfile.h" // // @@ -39,13 +38,17 @@ using namespace geom; Ewald::Ewald(StaticVals & stat, System & sys) : ff(stat.forcefield), mols(stat.mol), currentCoords(sys.coordinates), + currentCOM(sys.com), sysPotRef(sys.potential), lambdaRef(sys.lambdaRef), #ifdef VARIABLE_PARTICLE_NUMBER molLookup(sys.molLookup), #else molLookup(stat.molLookup), #endif - currentAxes(sys.boxDimRef), - currentCOM(sys.com), sysPotRef(sys.potential), lambdaRef(sys.lambdaRef) +#ifdef VARIABLE_VOLUME + currentAxes(sys.boxDimRef) +#else + currentAxes(*stat.GetBoxDim()) +#endif { ewald = false; electrostatic = false; @@ -100,6 +103,7 @@ Ewald::~Ewald() void Ewald::Init() { + for(uint m = 0; m < mols.count; ++m) { const MoleculeKind& molKind = mols.GetKind(m); for(uint a = 0; a < molKind.NumAtoms(); ++a) { @@ -118,13 +122,13 @@ void Ewald::Init() startMol.resize(currentCoords.Count()); lengthMol.resize(currentCoords.Count()); - for(int atom = 0; atom < (int) currentCoords.Count(); atom++) { + for(int atom = 0; atom < currentCoords.Count(); atom++) { startMol[atom] = mols.MolStart(particleMol[atom]); lengthMol[atom] = mols.MolLength(particleMol[atom]); } AllocMem(); - //initialize K vectors and reciprocal terms + //initialize K vectors and reciprocate terms UpdateVectorsAndRecipTerms(true); } @@ -191,20 +195,17 @@ void Ewald::AllocMem() } -//calculate reciprocal terms for a box. Should be called only at -//the start of the simulation to initialize the settings and when -//testing a change in box dimensions, such as a volume transfer. +//calculate reciprocal term for a box void Ewald::BoxReciprocalSetup(uint box, XYZArray const& molCoords) { if (box < BOXES_WITH_U_NB) { - GOMC_EVENT_START(1, GomcProfileEvent::RECIP_BOX_SETUP); MoleculeLookup::box_iterator end = molLookup.BoxEnd(box); MoleculeLookup::box_iterator thisMol = molLookup.BoxBegin(box); #ifdef GOMC_CUDA int numberOfAtoms = 0, i = 0; - for(int k = 0; k < (int) mols.GetKindsCount(); k++) { + for(int k = 0; k < mols.GetKindsCount(); k++) { MoleculeKind const& thisKind = mols.kinds[k]; numberOfAtoms += thisKind.NumAtoms() * molLookup.NumKindInBox(k, box); } @@ -242,10 +243,9 @@ void Ewald::BoxReciprocalSetup(uint box, XYZArray const& molCoords) uint start = mols.MolStart(*thisMol); #ifdef _OPENMP - #pragma omp parallel for default(none) shared(box, lambdaCoef, molCoords, \ - start, thisKind) + #pragma omp parallel for default(none) shared(box, lambdaCoef, molCoords, start, thisKind) #endif - for (int i = 0; i < (int) imageSize[box]; i++) { + for (int i = 0; i < imageSize[box]; i++) { double sumReal = 0.0; double sumImaginary = 0.0; @@ -270,133 +270,28 @@ void Ewald::BoxReciprocalSetup(uint box, XYZArray const& molCoords) thisMol++; } #endif - GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_BOX_SETUP); - } -} - - -//Calculate reciprocal terms for a box, when an updated value is needed -//because the number and location of molecules could have changed since -//the volume was set. Examples include MultiParticle and Molecule Exchange -//moves. For these calls, we need to use the Reference settings, since -//these hold the information for the current box dimensions. -void Ewald::BoxReciprocalSums(uint box, XYZArray const& molCoords) -{ - if (box < BOXES_WITH_U_NB) { - GOMC_EVENT_START(1, GomcProfileEvent::RECIP_BOX_SETUP); - MoleculeLookup::box_iterator thisMol = molLookup.BoxBegin(box); - MoleculeLookup::box_iterator end = molLookup.BoxEnd(box); - -#ifdef GOMC_CUDA - int numberOfAtoms = 0, i = 0; - - for(int k = 0; k < (int) mols.GetKindsCount(); k++) { - MoleculeKind const& thisKind = mols.kinds[k]; - numberOfAtoms += thisKind.NumAtoms() * molLookup.NumKindInBox(k, box); - } - - XYZArray thisBoxCoords(numberOfAtoms); - std::vector chargeBox; - - while (thisMol != end) { - MoleculeKind const& thisKind = mols.GetKind(*thisMol); - double lambdaCoef = GetLambdaCoef(*thisMol, box); - for (uint j = 0; j < thisKind.NumAtoms(); j++) { - thisBoxCoords.Set(i, molCoords[mols.MolStart(*thisMol) + j]); - chargeBox.push_back(thisKind.AtomCharge(j) * lambdaCoef); - i++; - } - thisMol++; - } - CallBoxReciprocalSumsGPU(ff.particles->getCUDAVars(), thisBoxCoords, - chargeBox, imageSizeRef[box], sumRnew[box], - sumInew[box], currentEnergyRecip[box], box); -#else -#ifdef _OPENMP - #pragma omp parallel default(none) shared(box) -#endif - { - std::memset(sumRnew[box], 0.0, sizeof(double) * imageSizeRef[box]); - std::memset(sumInew[box], 0.0, sizeof(double) * imageSizeRef[box]); - } - - while (thisMol != end) { - MoleculeKind const& thisKind = mols.GetKind(*thisMol); - double lambdaCoef = GetLambdaCoef(*thisMol, box); - uint startAtom = mols.MolStart(*thisMol); - -#ifdef _OPENMP - #pragma omp parallel for default(none) shared(box, lambdaCoef, molCoords, \ - startAtom, thisKind) -#endif - for (int i = 0; i < (int) imageSizeRef[box]; i++) { - double sumReal = 0.0; - double sumImaginary = 0.0; - - for (uint j = 0; j < thisKind.NumAtoms(); j++) { - unsigned long currentAtom = startAtom + j; - if(particleHasNoCharge[currentAtom]) { - continue; - } - double dotProduct = Dot(currentAtom, kxRef[box][i], kyRef[box][i], - kzRef[box][i], molCoords); - - // TODO: sincos() can be used to optimize (GNU compiler only) - // Windows doesn't have sincos() function and - // Intel compiler automatically optimizes this part - sumReal += (thisKind.AtomCharge(j) * cos(dotProduct)); - sumImaginary += (thisKind.AtomCharge(j) * sin(dotProduct)); - } - //we assume all atom charges are scaled with lambda - sumRnew[box][i] += (lambdaCoef * sumReal); - sumInew[box][i] += (lambdaCoef * sumImaginary); - } - thisMol++; - } -#endif - GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_BOX_SETUP); } } -//Calculate reciprocal energy for a box. This function is called for two reasons: -// 1. During a volume move, we call this function to total the sums for the new -// volume. For these calls, need to total the sums with the new volume settings. -// 2. During a Molecule Exchange or MultiParticle move, we need to recompute these -// sums for the current volume, since the number and location of molecules -// could have changed since the volume was set. For these calls, we need to -// use the Reference settings, since these hold the information for the current -// box dimensions. Also called at the start of the simulation, after the -// Reference volume parameters have been set. -double Ewald::BoxReciprocal(uint box, bool isNewVolume) const +//calculate reciprocal term for a box +double Ewald::BoxReciprocal(uint box) const { double energyRecip = 0.0; if (box < BOXES_WITH_U_NB) { - //GOMC_EVENT_START(1, GomcProfileEvent::RECIP_BOX_ENERGY); #ifdef GOMC_CUDA return currentEnergyRecip[box]; #else - double *prefactPtr; - int imageSzVal; - if (isNewVolume) { - prefactPtr = prefact[box]; - imageSzVal = static_cast(imageSize[box]); - } else { - prefactPtr = prefactRef[box]; - imageSzVal = static_cast(imageSizeRef[box]); - } #ifdef _OPENMP - #pragma omp parallel for default(none) shared(box, imageSzVal, prefactPtr) \ - reduction(+:energyRecip) + #pragma omp parallel for default(none) shared(box) reduction(+:energyRecip) #endif - for (int i = 0; i < imageSzVal; i++) { + for (int i = 0; i < imageSize[box]; i++) { energyRecip += (( sumRnew[box][i] * sumRnew[box][i] + sumInew[box][i] * sumInew[box][i]) * - prefactPtr[i]); + prefact[box][i]); } #endif - //GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_BOX_ENERGY); } return energyRecip; @@ -411,7 +306,6 @@ double Ewald::MolReciprocal(XYZArray const& molCoords, double energyRecipOld = 0.0; if (box < BOXES_WITH_U_NB) { - GOMC_EVENT_START(1, GomcProfileEvent::RECIP_MOL_ENERGY); MoleculeKind const& thisKind = mols.GetKind(molIndex); uint length = thisKind.NumAtoms(); uint startAtom = mols.MolStart(molIndex); @@ -438,7 +332,7 @@ reduction(+:energyRecipNew) reduction(+:energyRecipNew) #endif #endif - for (int i = 0; i < (int) imageSizeRef[box]; i++) { + for (int i = 0; i < imageSizeRef[box]; i++) { double sumRealNew = 0.0; double sumImaginaryNew = 0.0; double sumRealOld = 0.0; @@ -469,7 +363,6 @@ reduction(+:energyRecipNew) } #endif energyRecipOld = sysPotRef.boxEnergy[box].recip; - GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_MOL_ENERGY); } return energyRecipNew - energyRecipOld; } @@ -478,7 +371,7 @@ reduction(+:energyRecipNew) //calculate reciprocal term in destination box for swap move //No need to scale the charge with lambda, since this function will not be -// called in free energy of NeMTMC +// called in free energy of CFCMC double Ewald::SwapDestRecip(const cbmc::TrialMol &newMol, const uint box, const int molIndex) @@ -487,10 +380,10 @@ double Ewald::SwapDestRecip(const cbmc::TrialMol &newMol, double energyRecipOld = 0.0; if (box < BOXES_WITH_U_NB) { - GOMC_EVENT_START(1, GomcProfileEvent::RECIP_SWAP_ENERGY); MoleculeKind const& thisKind = newMol.GetKind(); XYZArray molCoords = newMol.GetCoords(); uint length = thisKind.NumAtoms(); + uint startAtom = mols.MolStart(molIndex); #ifdef GOMC_CUDA bool insert = true; std::vector MolCharge; @@ -502,7 +395,6 @@ double Ewald::SwapDestRecip(const cbmc::TrialMol &newMol, sumRnew[box], sumInew[box], insert, energyRecipNew, box); #else - uint startAtom = mols.MolStart(molIndex); #ifdef _OPENMP #if GCC_VERSION >= 90000 #pragma omp parallel for default(none) shared(length, molCoords, startAtom, \ @@ -512,7 +404,7 @@ thisKind, box) reduction(+:energyRecipNew) thisKind) reduction(+:energyRecipNew) #endif #endif - for (int i = 0; i < (int) imageSizeRef[box]; i++) { + for (int i = 0; i < imageSizeRef[box]; i++) { double sumRealNew = 0.0; double sumImaginaryNew = 0.0; @@ -536,36 +428,28 @@ thisKind) reduction(+:energyRecipNew) } #endif energyRecipOld = sysPotRef.boxEnergy[box].recip; - GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_SWAP_ENERGY); } return energyRecipNew - energyRecipOld; } //calculate reciprocal term for lambdaNew and Old with same coordinates -double Ewald::ChangeLambdaRecip(XYZArray const& molCoords, const double lambdaOld, - const double lambdaNew, const uint molIndex, - const uint box) +double Ewald::CFCMCRecip(XYZArray const& molCoords, const double lambdaOld, + const double lambdaNew, const uint molIndex, + const uint box) { double energyRecipNew = 0.0; double energyRecipOld = 0.0; + // + //Need to implement the GPU part + // if (box < BOXES_WITH_U_NB) { - GOMC_EVENT_START(1, GomcProfileEvent::RECIP_NEMTMC_ENERGY); MoleculeKind const& thisKind = mols.GetKind(molIndex); uint length = thisKind.NumAtoms(); uint startAtom = mols.MolStart(molIndex); double lambdaCoef = sqrt(lambdaNew) - sqrt(lambdaOld); -#ifdef GOMC_CUDA - std::vector MolCharge; - for(uint p = 0; p < length; p++) { - MolCharge.push_back(thisKind.AtomCharge(p)); - } - CallChangeLambdaMolReciprocalGPU(ff.particles->getCUDAVars(), - molCoords, MolCharge, imageSizeRef[box], - sumRnew[box], sumInew[box], energyRecipNew, - lambdaCoef, box); -#else + #ifdef _OPENMP #if GCC_VERSION >= 90000 #pragma omp parallel for default(none) shared(lambdaCoef, length, molCoords, \ @@ -577,7 +461,7 @@ reduction(+:energyRecipNew) reduction(+:energyRecipNew) #endif #endif - for (int i = 0; i < (int) imageSizeRef[box]; i++) { + for (int i = 0; i < imageSizeRef[box]; i++) { double sumRealNew = 0.0; double sumImaginaryNew = 0.0; @@ -601,9 +485,7 @@ reduction(+:energyRecipNew) energyRecipNew += (sumRnew[box][i] * sumRnew[box][i] + sumInew[box][i] * sumInew[box][i]) * prefactRef[box][i]; } -#endif energyRecipOld = sysPotRef.boxEnergy[box].recip; - GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_NEMTMC_ENERGY); } return energyRecipNew - energyRecipOld; @@ -634,7 +516,7 @@ reduction(+:energyRecip[:lambdaSize]) reduction(+:energyRecip[:lambdaSize]) #endif #endif - for (int i = 0; i < (int) imageSizeRef[box]; i++) { + for (uint i = 0; i < imageSizeRef[box]; i++) { double sumReal = 0.0; double sumImaginary = 0.0; @@ -671,19 +553,16 @@ reduction(+:energyRecip[:lambdaSize]) void Ewald::RecipInit(uint box, BoxDimensions const& boxAxes) { - GOMC_EVENT_START(1, GomcProfileEvent::RECIP_INIT); if(boxAxes.orthogonal[box]) RecipInitOrth(box, boxAxes); else RecipInitNonOrth(box, boxAxes); - - GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_INIT); } //calculate reciprocal term in source box for swap move //No need to scale the charge with lambda, since this function is not being -// called for free energy and NeMTMC +// called for free energy and CFCMC double Ewald::SwapSourceRecip(const cbmc::TrialMol &oldMol, const uint box, const int molIndex) { @@ -691,10 +570,10 @@ double Ewald::SwapSourceRecip(const cbmc::TrialMol &oldMol, double energyRecipOld = 0.0; if (box < BOXES_WITH_U_NB) { - GOMC_EVENT_START(1, GomcProfileEvent::RECIP_SWAP_ENERGY); MoleculeKind const& thisKind = oldMol.GetKind(); XYZArray molCoords = oldMol.GetCoords(); uint length = thisKind.NumAtoms(); + uint startAtom = mols.MolStart(molIndex); #ifdef GOMC_CUDA bool insert = false; std::vector MolCharge; @@ -707,7 +586,6 @@ double Ewald::SwapSourceRecip(const cbmc::TrialMol &oldMol, insert, energyRecipNew, box); #else - uint startAtom = mols.MolStart(molIndex); #ifdef _OPENMP #if GCC_VERSION >= 90000 #pragma omp parallel for default(none) shared(length, molCoords, startAtom, \ @@ -719,7 +597,7 @@ reduction(+:energyRecipNew) reduction(+:energyRecipNew) #endif #endif - for (int i = 0; i < (int) imageSizeRef[box]; i++) { + for (int i = 0; i < imageSizeRef[box]; i++) { double sumRealNew = 0.0; double sumImaginaryNew = 0.0; @@ -744,7 +622,6 @@ reduction(+:energyRecipNew) } #endif energyRecipOld = sysPotRef.boxEnergy[box].recip; - GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_SWAP_ENERGY); } return energyRecipNew - energyRecipOld; } @@ -752,11 +629,11 @@ reduction(+:energyRecipNew) //calculate reciprocal term for inserting some molecules (kindA) in destination // box and removing molecule (kindB) from destination box -double Ewald::MolExchangeReciprocal(const std::vector &newMol, - const std::vector &oldMol, - const std::vector &molIndexNew, - const std::vector &molIndexOld, - bool first_call) +double Ewald::SwapRecip(const std::vector &newMol, + const std::vector &oldMol, + const std::vector molIndexNew, + const std::vector molIndexOld, + bool first_call) { double energyRecipNew = 0.0; double energyRecipOld = 0.0; @@ -764,7 +641,6 @@ double Ewald::MolExchangeReciprocal(const std::vector &newMol, uint box = newMol[0].GetBox(); if (box < BOXES_WITH_U_NB) { - GOMC_EVENT_START(1, GomcProfileEvent::RECIP_MEMC_ENERGY); uint lengthNew, lengthOld; MoleculeKind const& thisKindNew = newMol[0].GetKind(); MoleculeKind const& thisKindOld = oldMol[0].GetKind(); @@ -772,11 +648,17 @@ double Ewald::MolExchangeReciprocal(const std::vector &newMol, lengthOld = thisKindOld.NumAtoms(); #ifdef _OPENMP +#if GCC_VERSION >= 90000 #pragma omp parallel for default(none) shared(box, first_call, lengthNew, lengthOld, \ newMol, oldMol, thisKindNew, thisKindOld, molIndexNew, molIndexOld) \ reduction(+:energyRecipNew) +#else + #pragma omp parallel for default(none) shared(box, first_call, lengthNew, lengthOld, \ + newMol, oldMol, thisKindNew, thisKindOld) \ +reduction(+:energyRecipNew) #endif - for (int i = 0; i < (int) imageSizeRef[box]; i++) { +#endif + for (int i = 0; i < imageSizeRef[box]; i++) { double sumRealNew = 0.0; double sumImaginaryNew = 0.0; @@ -847,19 +729,8 @@ reduction(+:energyRecipNew) // Keep hold of the old recip value energyRecipOld = sysPotRef.boxEnergy[box].recip; - GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_MEMC_ENERGY); } -//Because MolExchangeReciprocal does not have a matching GPU function, this is -//a stub function to copy the CPU sumRnew and sumInew vectors to the GPU in case -//the move is accepted. If this function is ported to the GPU, this call should -//be moved to the beginning of the MolExchangeReciprocal function and called -//instead of running the CPU code. -#ifdef GOMC_CUDA - CallMolExchangeReciprocalGPU(ff.particles->getCUDAVars(), imageSizeRef[box], - sumRnew[box], sumInew[box], box); -#endif - // Return the difference between old and new reciprocal energies return energyRecipNew - energyRecipOld; } @@ -901,12 +772,12 @@ void Ewald::RecipInitOrth(uint box, BoxDimensions const& boxAxes) double alpsqr4 = 1.0 / (4.0 * ff.alphaSq[box]); XYZ constValue = boxAxes.axis.Get(box); constValue.Inverse(); - constValue *= 2.0 * M_PI; + constValue *= 2 * M_PI; - double vol = boxAxes.volume[box] / (4.0 * M_PI); - nkx_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).x / (2.0 * M_PI)) + 1; - nky_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).y / (2.0 * M_PI)) + 1; - nkz_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).z / (2.0 * M_PI)) + 1; + double vol = boxAxes.volume[box] / (4 * M_PI); + nkx_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).x / (2 * M_PI)) + 1; + nky_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).y / (2 * M_PI)) + 1; + nkz_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).z / (2 * M_PI)) + 1; kmax[box] = std::max(std::max(nkx_max, nky_max), std::max(nky_max, nkz_max)); for(x = 0; x <= nkx_max; x++) { @@ -961,12 +832,12 @@ void Ewald::RecipInitNonOrth(uint box, BoxDimensions const& boxAxes) cellB.Scale(2, boxAxes.axis.Get(box).z); XYZArray cellB_Inv(3); double det = cellB.AdjointMatrix(cellB_Inv); - cellB_Inv.ScaleRange(0, 3, (2.0 * M_PI) / det); + cellB_Inv.ScaleRange(0, 3, (2 * M_PI) / det); - double vol = boxAxes.volume[box] / (4.0 * M_PI); - nkx_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).x / (2.0 * M_PI)) + 1; - nky_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).y / (2.0 * M_PI)) + 1; - nkz_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).z / (2.0 * M_PI)) + 1; + double vol = boxAxes.volume[box] / (4 * M_PI); + nkx_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).x / (2 * M_PI)) + 1; + nky_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).y / (2 * M_PI)) + 1; + nkz_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).z / (2 * M_PI)) + 1; kmax[box] = std::max(std::max(nkx_max, nky_max), std::max(nky_max, nkz_max)); for (x = 0; x <= nkx_max; x++) { @@ -1030,11 +901,11 @@ void Ewald::RecipCountInit(uint box, BoxDimensions const& boxAxes) cellB.Scale(2, constValue.z); XYZArray cellB_Inv(3); double det = cellB.AdjointMatrix(cellB_Inv); - cellB_Inv.ScaleRange(0, 3, (2.0 * M_PI) / det); + cellB_Inv.ScaleRange(0, 3, (2 * M_PI) / det); - nkx_max = int(ff.recip_rcut[box] * constValue.x / (2.0 * M_PI)) + 1; - nky_max = int(ff.recip_rcut[box] * constValue.y / (2.0 * M_PI)) + 1; - nkz_max = int(ff.recip_rcut[box] * constValue.z / (2.0 * M_PI)) + 1; + nkx_max = int(ff.recip_rcut[box] * constValue.x / (2 * M_PI)) + 1; + nky_max = int(ff.recip_rcut[box] * constValue.y / (2 * M_PI)) + 1; + nkz_max = int(ff.recip_rcut[box] * constValue.z / (2 * M_PI)) + 1; for(x = 0; x <= nkx_max; x++) { if(x == 0.0) @@ -1079,9 +950,12 @@ void Ewald::SetRecipRef(uint box) std::memcpy(prefactRef[box], prefact[box], sizeof(double) *imageSize[box]); } #ifdef GOMC_CUDA - CopyCurrentToRefCUDA(ff.particles->getCUDAVars(), box, imageSize[box]); + CopyCurrentToRefCUDA(ff.particles->getCUDAVars(), + box, imageSize[box]); #endif - imageSizeRef[box] = imageSize[box]; + for(uint b = 0; b < BOXES_WITH_U_NB; b++) { + imageSizeRef[b] = imageSize[b]; + } } //calculate correction term for a molecule, with system lambda @@ -1090,7 +964,6 @@ double Ewald::MolCorrection(uint molIndex, uint box) const if (box >= BOXES_WITH_U_NB) return 0.0; - GOMC_EVENT_START(1, GomcProfileEvent::CORR_MOL); double dist, distSq; double correction = 0.0; XYZ virComponents; @@ -1105,21 +978,14 @@ double Ewald::MolCorrection(uint molIndex, uint box) const continue; } for (uint j = i + 1; j < atomSize; j++) { - if(currentAxes.InRcut(distSq, virComponents, currentCoords, + currentAxes.InRcut(distSq, virComponents, currentCoords, start + i, start + j, box); dist = sqrt(distSq); correction += (thisKind.AtomCharge(i) * thisKind.AtomCharge(j) * erf(ff.alpha[box] * dist) / dist); } - // within cutoff - // CalcCoul = qiqj/r(erfc(alpha*r)+erf(alpha*r)) - - // outside cutoff - // erf(alpha*r)*qiqj/r - // CalcCoul = erf } - GOMC_EVENT_STOP(1, GomcProfileEvent::CORR_MOL); return -1.0 * num::qqFact * correction * lambdaCoef * lambdaCoef; } @@ -1162,12 +1028,11 @@ void Ewald::ChangeCorrection(Energy *energyDiff, Energy &dUdL_Coul, } //calculate self term for a box, using system lambda -double Ewald::BoxSelf(uint box) const +double Ewald::BoxSelf(BoxDimensions const& boxAxes, uint box) const { if (box >= BOXES_WITH_U_NB) return 0.0; - GOMC_EVENT_START(1, GomcProfileEvent::SELF_BOX); double self = 0.0; double molSelfEnergy; uint i, j, length, molNum; @@ -1196,10 +1061,8 @@ double Ewald::BoxSelf(uint box) const } } - // M_2_SQRTPI is 2/sqrt(PI), so need to multiply by 0.5 to get sqrt(PI) - self *= -1.0 * ff.alpha[box] * num::qqFact * M_2_SQRTPI * 0.5; + self = -1.0 * self * ff.alpha[box] * num::qqFact / sqrt(M_PI); - GOMC_EVENT_STOP(1, GomcProfileEvent::SELF_BOX); return self; } @@ -1212,12 +1075,12 @@ Virial Ewald::VirialReciprocal(Virial& virial, uint box) const if (box >= BOXES_WITH_U_NB) return tempVir; - GOMC_EVENT_START(1, GomcProfileEvent::RECIP_BOX_VIRIAL); double wT11 = 0.0, wT12 = 0.0, wT13 = 0.0; double wT22 = 0.0, wT23 = 0.0, wT33 = 0.0; + double recipIntra = 0.0; double constVal = 1.0 / (4.0 * ff.alphaSq[box]); - double lambdaCoef; + double charge, lambdaCoef; uint p, length, startAtom, atom; MoleculeLookup::box_iterator thisMol = molLookup.BoxBegin(box), @@ -1228,7 +1091,7 @@ Virial Ewald::VirialReciprocal(Virial& virial, uint box) const #ifdef GOMC_CUDA int numberOfAtoms = 0, atomIndex = 0; - for(int k = 0; k < (int) mols.GetKindsCount(); k++) { + for(int k = 0; k < mols.GetKindsCount(); k++) { MoleculeKind const& thisKind = mols.kinds[k]; numberOfAtoms += thisKind.NumAtoms() * molLookup.NumKindInBox(k, box); } @@ -1268,7 +1131,7 @@ Virial Ewald::VirialReciprocal(Virial& virial, uint box) const #ifdef _OPENMP #pragma omp parallel for default(none) shared(box, constVal) reduction(+:wT11, wT22, wT33) #endif - for (int i = 0; i < (int) imageSizeRef[box]; i++) { + for (int i = 0; i < imageSizeRef[box]; i++) { double factor = prefactRef[box][i] * (sumRref[box][i] * sumRref[box][i] + sumIref[box][i] * sumIref[box][i]); @@ -1300,12 +1163,12 @@ Virial Ewald::VirialReciprocal(Virial& virial, uint box) const currentAxes.UnwrapPBC(atomC, box, comC); diffC = atomC - comC; //scale the charge with lambda for Free energy calc - double charge = particleCharge[atom] * lambdaCoef; + charge = particleCharge[atom] * lambdaCoef; #ifdef _OPENMP #pragma omp parallel for default(none) shared(atom, box, charge, diffC) reduction(+:wT11, wT22, wT33) #endif - for (int i = 0; i < (int) imageSizeRef[box]; i++) { + for (int i = 0; i < imageSizeRef[box]; i++) { //compute the dot product of k and r double arg = Dot(atom, kxRef[box][i], kyRef[box][i], kzRef[box][i], currentCoords); @@ -1339,7 +1202,6 @@ Virial Ewald::VirialReciprocal(Virial& virial, uint box) const // setting virial of reciprocal space tempVir.recip = wT11 + wT22 + wT33; - GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_BOX_VIRIAL); return tempVir; } @@ -1354,7 +1216,6 @@ double Ewald::SwapCorrection(const cbmc::TrialMol& trialMol) const if (box >= BOXES_WITH_U_NB) return 0.0; - GOMC_EVENT_START(1, GomcProfileEvent::CORR_SWAP); double dist, distSq; double correction = 0.0; XYZ virComponents; @@ -1371,7 +1232,6 @@ double Ewald::SwapCorrection(const cbmc::TrialMol& trialMol) const erf(ff.alpha[box] * dist) / dist); } } - GOMC_EVENT_STOP(1, GomcProfileEvent::CORR_SWAP); return num::qqFact * correction; } @@ -1384,7 +1244,6 @@ double Ewald::SwapCorrection(const cbmc::TrialMol& trialMol, if (box >= BOXES_WITH_U_NB) return 0.0; - GOMC_EVENT_START(1, GomcProfileEvent::CORR_SWAP); double dist, distSq; double correction = 0.0; XYZ virComponents; @@ -1406,20 +1265,18 @@ double Ewald::SwapCorrection(const cbmc::TrialMol& trialMol, erf(ff.alpha[box] * dist) / dist); } } - GOMC_EVENT_STOP(1, GomcProfileEvent::CORR_SWAP); return num::qqFact * correction * lambdaCoef * lambdaCoef; } //It's called if we transfer one molecule from one box to another //No need to scale the charge with lambda, since this function is not being -// called from free energy or NeMTMC +// called from free energy or CFCMC double Ewald::SwapSelf(const cbmc::TrialMol& trialMol) const { uint box = trialMol.GetBox(); if (box >= BOXES_WITH_U_NB) return 0.0; - GOMC_EVENT_START(1, GomcProfileEvent::SELF_SWAP); MoleculeKind const& thisKind = trialMol.GetKind(); uint atomSize = thisKind.NumAtoms(); double en_self = 0.0; @@ -1427,9 +1284,7 @@ double Ewald::SwapSelf(const cbmc::TrialMol& trialMol) const for (uint i = 0; i < atomSize; i++) { en_self -= (thisKind.AtomCharge(i) * thisKind.AtomCharge(i)); } - GOMC_EVENT_STOP(1, GomcProfileEvent::SELF_SWAP); - // M_2_SQRTPI is 2/sqrt(PI), so need to multiply by 0.5 to get sqrt(PI) - return (en_self * ff.alpha[box] * num::qqFact * M_2_SQRTPI * 0.5); + return (en_self * ff.alpha[box] * num::qqFact / sqrt(M_PI)); } //It's called in free energy calculation to calculate the change in @@ -1447,8 +1302,7 @@ void Ewald::ChangeSelf(Energy *energyDiff, Energy &dUdL_Coul, for (uint i = 0; i < atomSize; i++) { en_self += (particleCharge[i + start] * particleCharge[i + start]); } - // M_2_SQRTPI is 2/sqrt(PI), so need to multiply by 0.5 to get sqrt(PI) - en_self *= -1.0 * ff.alpha[box] * num::qqFact * M_2_SQRTPI * 0.5; + en_self *= -1.0 * ff.alpha[box] * num::qqFact / sqrt(M_PI); //Calculate the energy difference for each lambda state for (uint s = 0; s < lambdaSize; s++) { @@ -1477,28 +1331,6 @@ void Ewald::UpdateRecip(uint box) #endif } -//copy reciprocal values from ref to new -//Used to reinitialize these two arrays for MultiParticle moves so that the -//forces can be computed using the current positions, before making the MP -//move. The other option would be to create different versions of these -//functions to access SumIref and SumRref instead of SumInew and SumRnew. -void Ewald::CopyRecip(uint box) -{ - if (box >= BOXES_WITH_U_NB) - return; - -#ifdef _OPENMP - #pragma omp parallel default(none) shared(box) -#endif - { - std::memcpy(sumRnew[box], sumRref[box], sizeof(double) * imageSizeRef[box]); - std::memcpy(sumInew[box], sumIref[box], sizeof(double) * imageSizeRef[box]); - } -#ifdef GOMC_CUDA - CopyRefToNewCUDA(ff.particles->getCUDAVars(), box, imageSizeRef[box]); -#endif -} - void Ewald::UpdateRecipVec(uint box) { double *tempKx, *tempKy, *tempKz, *tempHsqr, *tempPrefact; @@ -1519,12 +1351,13 @@ void Ewald::UpdateRecipVec(uint box) kz[box] = tempKz; hsqr[box] = tempHsqr; prefact[box] = tempPrefact; - - imageSizeRef[box] = imageSize[box]; - #ifdef GOMC_CUDA UpdateRecipVecCUDA(ff.particles->getCUDAVars(), box); #endif + + for(uint b = 0; b < BOXES_WITH_U_NB; b++) { + imageSizeRef[b] = imageSize[b]; + } } void compareDouble(const double &x, const double &y, const int &i) @@ -1542,35 +1375,23 @@ void Ewald::BoxForceReciprocal(XYZArray const& molCoords, uint box) { if(multiParticleEnabled && (box < BOXES_WITH_U_NB)) { - GOMC_EVENT_START(1, GomcProfileEvent::RECIP_BOX_FORCE); - // M_2_SQRTPI is 2/sqrt(PI) - double constValue = ff.alpha[box] * M_2_SQRTPI; + double constValue = 2.0 * ff.alpha[box] / sqrt(M_PI); #ifdef GOMC_CUDA - bool *particleUsed; - particleUsed = new bool[atomForceRec.Count()]; - memset((void *) particleUsed, false, atomForceRec.Count() * sizeof(bool)); -#if ENSEMBLE == GEMC || ENSEMBLE == GCMC - memset((void *) particleUsed, false, atomForceRec.Count() * sizeof(bool)); MoleculeLookup::box_iterator thisMol = molLookup.BoxBegin(box); MoleculeLookup::box_iterator end = molLookup.BoxEnd(box); while(thisMol != end) { uint molIndex = *thisMol; uint length = mols.GetKind(molIndex).NumAtoms(); uint start = mols.MolStart(molIndex); - for(uint p = 0; p < length; p++) { + for(int p = 0; p < length; p++) { atomForceRec.Set(start + p, 0.0, 0.0, 0.0); - particleUsed[start + p] = true; } molForceRec.Set(molIndex, 0.0, 0.0, 0.0); thisMol++; } -#else - //Only one box, so clear all atoms and molecules and mark all particles as Used - atomForceRec.Reset(); - molForceRec.Reset(); - memset((void *) particleUsed, true, atomForceRec.Count() * sizeof(bool)); -#endif + // initialize the start and end of each box + initializeBoxRange(); CallBoxForceReciprocalGPU( ff.particles->getCUDAVars(), @@ -1578,20 +1399,21 @@ void Ewald::BoxForceReciprocal(XYZArray const& molCoords, molForceRec, particleCharge, particleMol, + particleKind, particleHasNoCharge, - particleUsed, startMol, lengthMol, ff.alpha[box], ff.alphaSq[box], num::qqFact, constValue, - imageSizeRef[box], + imageSize[box], molCoords, + boxStart[box], + boxEnd[box], currentAxes, box ); - delete[] particleUsed; #else // molecule iterator MoleculeLookup::box_iterator thisMol = molLookup.BoxBegin(box); @@ -1609,6 +1431,7 @@ void Ewald::BoxForceReciprocal(XYZArray const& molCoords, for(p = start; p < start + length; p++) { double X = 0.0, Y = 0.0, Z = 0.0; + double intraForce; if(!particleHasNoCharge[p]) { // subtract the intra forces(correction) @@ -1619,7 +1442,7 @@ void Ewald::BoxForceReciprocal(XYZArray const& molCoords, double dist = sqrt(distSq); double expConstValue = exp(-1.0 * ff.alphaSq[box] * distSq); double qiqj = particleCharge[p] * particleCharge[j] * num::qqFact; - double intraForce = qiqj * lambdaCoef * lambdaCoef / distSq; + intraForce = qiqj * lambdaCoef * lambdaCoef / distSq; intraForce *= ((erf(ff.alpha[box] * dist) / dist) - constValue * expConstValue); X -= intraForce * distVect.x; @@ -1630,15 +1453,15 @@ void Ewald::BoxForceReciprocal(XYZArray const& molCoords, #ifdef _OPENMP #pragma omp parallel for default(none) shared(box, lambdaCoef, molCoords, p) reduction(+:X, Y, Z) #endif - for(int i = 0; i < (int) imageSizeRef[box]; i++) { - double dot = Dot(p, kxRef[box][i], kyRef[box][i], kzRef[box][i], molCoords); + for(int i = 0; i < imageSize[box]; i++) { + double dot = Dot(p, kx[box][i], ky[box][i], kz[box][i], molCoords); - double factor = 2.0 * particleCharge[p] * prefactRef[box][i] * lambdaCoef * + double factor = 2.0 * particleCharge[p] * prefact[box][i] * lambdaCoef * (sin(dot) * sumRnew[box][i] - cos(dot) * sumInew[box][i]); - X += factor * kxRef[box][i]; - Y += factor * kyRef[box][i]; - Z += factor * kzRef[box][i]; + X += factor * kx[box][i]; + Y += factor * ky[box][i]; + Z += factor * kz[box][i]; } } atomForceRec.Set(p, X, Y, Z); @@ -1647,13 +1470,28 @@ void Ewald::BoxForceReciprocal(XYZArray const& molCoords, thisMol++; } #endif - GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_BOX_FORCE); } } double Ewald::GetLambdaCoef(uint molA, uint box) const { - double lambda = lambdaRef.GetLambdaCoulomb(molA, box); + double lambda = lambdaRef.GetLambdaCoulomb(molA, mols.GetMolKind(molA), box); //Each charge gets sq root of it. return sqrt(lambda); } + +void Ewald::initializeBoxRange() +{ + if(BOX_TOTAL == 1) { + boxStart[mv::BOX0] = 0; + boxEnd[mv::BOX0] = currentCoords.Count(); + } else { + boxStart[mv::BOX0] = 0; + MoleculeLookup::box_iterator startBox1 = molLookup.BoxBegin(mv::BOX1); + int firstMoleculeInBox1 = *startBox1; + int indexOfFirstAtomInBox1 = mols.MolStart(firstMoleculeInBox1); + boxEnd[mv::BOX0] = indexOfFirstAtomInBox1; + boxStart[mv::BOX1] = indexOfFirstAtomInBox1; + boxEnd[mv::BOX1] = currentCoords.Count(); + } +} From 074d28167610d19a56fac0898be9c0e4eb9ceb52 Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Mon, 22 Nov 2021 19:03:58 -0600 Subject: [PATCH 06/16] Revert "change ewald back" This reverts commit 047dc3cf40d6b8684ee07d1933707de32ccea6e2. --- src/Ewald.cpp | 390 +++++++++++++++++++++++++++++++++++--------------- 1 file changed, 276 insertions(+), 114 deletions(-) diff --git a/src/Ewald.cpp b/src/Ewald.cpp index 6895148ed..8d19fa151 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -24,6 +24,7 @@ along with this program, also can be found at . #include "CalculateForceCUDAKernel.cuh" #include "ConstantDefinitionsCUDAKernel.cuh" #endif +#include "GOMCEventsProfile.h" // // @@ -38,17 +39,13 @@ using namespace geom; Ewald::Ewald(StaticVals & stat, System & sys) : ff(stat.forcefield), mols(stat.mol), currentCoords(sys.coordinates), - currentCOM(sys.com), sysPotRef(sys.potential), lambdaRef(sys.lambdaRef), #ifdef VARIABLE_PARTICLE_NUMBER molLookup(sys.molLookup), #else molLookup(stat.molLookup), #endif -#ifdef VARIABLE_VOLUME - currentAxes(sys.boxDimRef) -#else - currentAxes(*stat.GetBoxDim()) -#endif + currentAxes(sys.boxDimRef), + currentCOM(sys.com), sysPotRef(sys.potential), lambdaRef(sys.lambdaRef) { ewald = false; electrostatic = false; @@ -103,7 +100,6 @@ Ewald::~Ewald() void Ewald::Init() { - for(uint m = 0; m < mols.count; ++m) { const MoleculeKind& molKind = mols.GetKind(m); for(uint a = 0; a < molKind.NumAtoms(); ++a) { @@ -122,13 +118,13 @@ void Ewald::Init() startMol.resize(currentCoords.Count()); lengthMol.resize(currentCoords.Count()); - for(int atom = 0; atom < currentCoords.Count(); atom++) { + for(int atom = 0; atom < (int) currentCoords.Count(); atom++) { startMol[atom] = mols.MolStart(particleMol[atom]); lengthMol[atom] = mols.MolLength(particleMol[atom]); } AllocMem(); - //initialize K vectors and reciprocate terms + //initialize K vectors and reciprocal terms UpdateVectorsAndRecipTerms(true); } @@ -195,17 +191,20 @@ void Ewald::AllocMem() } -//calculate reciprocal term for a box +//calculate reciprocal terms for a box. Should be called only at +//the start of the simulation to initialize the settings and when +//testing a change in box dimensions, such as a volume transfer. void Ewald::BoxReciprocalSetup(uint box, XYZArray const& molCoords) { if (box < BOXES_WITH_U_NB) { + GOMC_EVENT_START(1, GomcProfileEvent::RECIP_BOX_SETUP); MoleculeLookup::box_iterator end = molLookup.BoxEnd(box); MoleculeLookup::box_iterator thisMol = molLookup.BoxBegin(box); #ifdef GOMC_CUDA int numberOfAtoms = 0, i = 0; - for(int k = 0; k < mols.GetKindsCount(); k++) { + for(int k = 0; k < (int) mols.GetKindsCount(); k++) { MoleculeKind const& thisKind = mols.kinds[k]; numberOfAtoms += thisKind.NumAtoms() * molLookup.NumKindInBox(k, box); } @@ -243,9 +242,10 @@ void Ewald::BoxReciprocalSetup(uint box, XYZArray const& molCoords) uint start = mols.MolStart(*thisMol); #ifdef _OPENMP - #pragma omp parallel for default(none) shared(box, lambdaCoef, molCoords, start, thisKind) + #pragma omp parallel for default(none) shared(box, lambdaCoef, molCoords, \ + start, thisKind) #endif - for (int i = 0; i < imageSize[box]; i++) { + for (int i = 0; i < (int) imageSize[box]; i++) { double sumReal = 0.0; double sumImaginary = 0.0; @@ -270,28 +270,133 @@ void Ewald::BoxReciprocalSetup(uint box, XYZArray const& molCoords) thisMol++; } #endif + GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_BOX_SETUP); + } +} + + +//Calculate reciprocal terms for a box, when an updated value is needed +//because the number and location of molecules could have changed since +//the volume was set. Examples include MultiParticle and Molecule Exchange +//moves. For these calls, we need to use the Reference settings, since +//these hold the information for the current box dimensions. +void Ewald::BoxReciprocalSums(uint box, XYZArray const& molCoords) +{ + if (box < BOXES_WITH_U_NB) { + GOMC_EVENT_START(1, GomcProfileEvent::RECIP_BOX_SETUP); + MoleculeLookup::box_iterator thisMol = molLookup.BoxBegin(box); + MoleculeLookup::box_iterator end = molLookup.BoxEnd(box); + +#ifdef GOMC_CUDA + int numberOfAtoms = 0, i = 0; + + for(int k = 0; k < (int) mols.GetKindsCount(); k++) { + MoleculeKind const& thisKind = mols.kinds[k]; + numberOfAtoms += thisKind.NumAtoms() * molLookup.NumKindInBox(k, box); + } + + XYZArray thisBoxCoords(numberOfAtoms); + std::vector chargeBox; + + while (thisMol != end) { + MoleculeKind const& thisKind = mols.GetKind(*thisMol); + double lambdaCoef = GetLambdaCoef(*thisMol, box); + for (uint j = 0; j < thisKind.NumAtoms(); j++) { + thisBoxCoords.Set(i, molCoords[mols.MolStart(*thisMol) + j]); + chargeBox.push_back(thisKind.AtomCharge(j) * lambdaCoef); + i++; + } + thisMol++; + } + CallBoxReciprocalSumsGPU(ff.particles->getCUDAVars(), thisBoxCoords, + chargeBox, imageSizeRef[box], sumRnew[box], + sumInew[box], currentEnergyRecip[box], box); +#else +#ifdef _OPENMP + #pragma omp parallel default(none) shared(box) +#endif + { + std::memset(sumRnew[box], 0.0, sizeof(double) * imageSizeRef[box]); + std::memset(sumInew[box], 0.0, sizeof(double) * imageSizeRef[box]); + } + + while (thisMol != end) { + MoleculeKind const& thisKind = mols.GetKind(*thisMol); + double lambdaCoef = GetLambdaCoef(*thisMol, box); + uint startAtom = mols.MolStart(*thisMol); + +#ifdef _OPENMP + #pragma omp parallel for default(none) shared(box, lambdaCoef, molCoords, \ + startAtom, thisKind) +#endif + for (int i = 0; i < (int) imageSizeRef[box]; i++) { + double sumReal = 0.0; + double sumImaginary = 0.0; + + for (uint j = 0; j < thisKind.NumAtoms(); j++) { + unsigned long currentAtom = startAtom + j; + if(particleHasNoCharge[currentAtom]) { + continue; + } + double dotProduct = Dot(currentAtom, kxRef[box][i], kyRef[box][i], + kzRef[box][i], molCoords); + + // TODO: sincos() can be used to optimize (GNU compiler only) + // Windows doesn't have sincos() function and + // Intel compiler automatically optimizes this part + sumReal += (thisKind.AtomCharge(j) * cos(dotProduct)); + sumImaginary += (thisKind.AtomCharge(j) * sin(dotProduct)); + } + //we assume all atom charges are scaled with lambda + sumRnew[box][i] += (lambdaCoef * sumReal); + sumInew[box][i] += (lambdaCoef * sumImaginary); + } + thisMol++; + } +#endif + GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_BOX_SETUP); } } -//calculate reciprocal term for a box -double Ewald::BoxReciprocal(uint box) const +//Calculate reciprocal energy for a box. This function is called for two reasons: +// 1. During a volume move, we call this function to total the sums for the new +// volume. For these calls, need to total the sums with the new volume settings. +// 2. During a Molecule Exchange or MultiParticle move, we need to recompute these +// sums for the current volume, since the number and location of molecules +// could have changed since the volume was set. For these calls, we need to +// use the Reference settings, since these hold the information for the current +// box dimensions. Also called at the start of the simulation, after the +// Reference volume parameters have been set. +double Ewald::BoxReciprocal(uint box, bool isNewVolume) const { double energyRecip = 0.0; if (box < BOXES_WITH_U_NB) { + //GOMC_EVENT_START(1, GomcProfileEvent::RECIP_BOX_ENERGY); #ifdef GOMC_CUDA return currentEnergyRecip[box]; #else + double *prefactPtr; + int imageSzVal; + if (isNewVolume) { + prefactPtr = prefact[box]; + imageSzVal = static_cast(imageSize[box]); + } else { + prefactPtr = prefactRef[box]; + imageSzVal = static_cast(imageSizeRef[box]); + } #ifdef _OPENMP - #pragma omp parallel for default(none) shared(box) reduction(+:energyRecip) + #pragma omp parallel for default(none) shared(box, imageSzVal, prefactPtr) \ + reduction(+:energyRecip) #endif - for (int i = 0; i < imageSize[box]; i++) { + for (int i = 0; i < imageSzVal; i++) { energyRecip += (( sumRnew[box][i] * sumRnew[box][i] + sumInew[box][i] * sumInew[box][i]) * - prefact[box][i]); + prefactPtr[i]); } #endif + //GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_BOX_ENERGY); } return energyRecip; @@ -306,6 +411,7 @@ double Ewald::MolReciprocal(XYZArray const& molCoords, double energyRecipOld = 0.0; if (box < BOXES_WITH_U_NB) { + GOMC_EVENT_START(1, GomcProfileEvent::RECIP_MOL_ENERGY); MoleculeKind const& thisKind = mols.GetKind(molIndex); uint length = thisKind.NumAtoms(); uint startAtom = mols.MolStart(molIndex); @@ -332,7 +438,7 @@ reduction(+:energyRecipNew) reduction(+:energyRecipNew) #endif #endif - for (int i = 0; i < imageSizeRef[box]; i++) { + for (int i = 0; i < (int) imageSizeRef[box]; i++) { double sumRealNew = 0.0; double sumImaginaryNew = 0.0; double sumRealOld = 0.0; @@ -363,6 +469,7 @@ reduction(+:energyRecipNew) } #endif energyRecipOld = sysPotRef.boxEnergy[box].recip; + GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_MOL_ENERGY); } return energyRecipNew - energyRecipOld; } @@ -371,7 +478,7 @@ reduction(+:energyRecipNew) //calculate reciprocal term in destination box for swap move //No need to scale the charge with lambda, since this function will not be -// called in free energy of CFCMC +// called in free energy of NeMTMC double Ewald::SwapDestRecip(const cbmc::TrialMol &newMol, const uint box, const int molIndex) @@ -380,10 +487,10 @@ double Ewald::SwapDestRecip(const cbmc::TrialMol &newMol, double energyRecipOld = 0.0; if (box < BOXES_WITH_U_NB) { + GOMC_EVENT_START(1, GomcProfileEvent::RECIP_SWAP_ENERGY); MoleculeKind const& thisKind = newMol.GetKind(); XYZArray molCoords = newMol.GetCoords(); uint length = thisKind.NumAtoms(); - uint startAtom = mols.MolStart(molIndex); #ifdef GOMC_CUDA bool insert = true; std::vector MolCharge; @@ -395,6 +502,7 @@ double Ewald::SwapDestRecip(const cbmc::TrialMol &newMol, sumRnew[box], sumInew[box], insert, energyRecipNew, box); #else + uint startAtom = mols.MolStart(molIndex); #ifdef _OPENMP #if GCC_VERSION >= 90000 #pragma omp parallel for default(none) shared(length, molCoords, startAtom, \ @@ -404,7 +512,7 @@ thisKind, box) reduction(+:energyRecipNew) thisKind) reduction(+:energyRecipNew) #endif #endif - for (int i = 0; i < imageSizeRef[box]; i++) { + for (int i = 0; i < (int) imageSizeRef[box]; i++) { double sumRealNew = 0.0; double sumImaginaryNew = 0.0; @@ -428,28 +536,36 @@ thisKind) reduction(+:energyRecipNew) } #endif energyRecipOld = sysPotRef.boxEnergy[box].recip; + GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_SWAP_ENERGY); } return energyRecipNew - energyRecipOld; } //calculate reciprocal term for lambdaNew and Old with same coordinates -double Ewald::CFCMCRecip(XYZArray const& molCoords, const double lambdaOld, - const double lambdaNew, const uint molIndex, - const uint box) +double Ewald::ChangeLambdaRecip(XYZArray const& molCoords, const double lambdaOld, + const double lambdaNew, const uint molIndex, + const uint box) { double energyRecipNew = 0.0; double energyRecipOld = 0.0; - // - //Need to implement the GPU part - // if (box < BOXES_WITH_U_NB) { + GOMC_EVENT_START(1, GomcProfileEvent::RECIP_NEMTMC_ENERGY); MoleculeKind const& thisKind = mols.GetKind(molIndex); uint length = thisKind.NumAtoms(); uint startAtom = mols.MolStart(molIndex); double lambdaCoef = sqrt(lambdaNew) - sqrt(lambdaOld); - +#ifdef GOMC_CUDA + std::vector MolCharge; + for(uint p = 0; p < length; p++) { + MolCharge.push_back(thisKind.AtomCharge(p)); + } + CallChangeLambdaMolReciprocalGPU(ff.particles->getCUDAVars(), + molCoords, MolCharge, imageSizeRef[box], + sumRnew[box], sumInew[box], energyRecipNew, + lambdaCoef, box); +#else #ifdef _OPENMP #if GCC_VERSION >= 90000 #pragma omp parallel for default(none) shared(lambdaCoef, length, molCoords, \ @@ -461,7 +577,7 @@ reduction(+:energyRecipNew) reduction(+:energyRecipNew) #endif #endif - for (int i = 0; i < imageSizeRef[box]; i++) { + for (int i = 0; i < (int) imageSizeRef[box]; i++) { double sumRealNew = 0.0; double sumImaginaryNew = 0.0; @@ -485,7 +601,9 @@ reduction(+:energyRecipNew) energyRecipNew += (sumRnew[box][i] * sumRnew[box][i] + sumInew[box][i] * sumInew[box][i]) * prefactRef[box][i]; } +#endif energyRecipOld = sysPotRef.boxEnergy[box].recip; + GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_NEMTMC_ENERGY); } return energyRecipNew - energyRecipOld; @@ -516,7 +634,7 @@ reduction(+:energyRecip[:lambdaSize]) reduction(+:energyRecip[:lambdaSize]) #endif #endif - for (uint i = 0; i < imageSizeRef[box]; i++) { + for (int i = 0; i < (int) imageSizeRef[box]; i++) { double sumReal = 0.0; double sumImaginary = 0.0; @@ -553,16 +671,19 @@ reduction(+:energyRecip[:lambdaSize]) void Ewald::RecipInit(uint box, BoxDimensions const& boxAxes) { + GOMC_EVENT_START(1, GomcProfileEvent::RECIP_INIT); if(boxAxes.orthogonal[box]) RecipInitOrth(box, boxAxes); else RecipInitNonOrth(box, boxAxes); + + GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_INIT); } //calculate reciprocal term in source box for swap move //No need to scale the charge with lambda, since this function is not being -// called for free energy and CFCMC +// called for free energy and NeMTMC double Ewald::SwapSourceRecip(const cbmc::TrialMol &oldMol, const uint box, const int molIndex) { @@ -570,10 +691,10 @@ double Ewald::SwapSourceRecip(const cbmc::TrialMol &oldMol, double energyRecipOld = 0.0; if (box < BOXES_WITH_U_NB) { + GOMC_EVENT_START(1, GomcProfileEvent::RECIP_SWAP_ENERGY); MoleculeKind const& thisKind = oldMol.GetKind(); XYZArray molCoords = oldMol.GetCoords(); uint length = thisKind.NumAtoms(); - uint startAtom = mols.MolStart(molIndex); #ifdef GOMC_CUDA bool insert = false; std::vector MolCharge; @@ -586,6 +707,7 @@ double Ewald::SwapSourceRecip(const cbmc::TrialMol &oldMol, insert, energyRecipNew, box); #else + uint startAtom = mols.MolStart(molIndex); #ifdef _OPENMP #if GCC_VERSION >= 90000 #pragma omp parallel for default(none) shared(length, molCoords, startAtom, \ @@ -597,7 +719,7 @@ reduction(+:energyRecipNew) reduction(+:energyRecipNew) #endif #endif - for (int i = 0; i < imageSizeRef[box]; i++) { + for (int i = 0; i < (int) imageSizeRef[box]; i++) { double sumRealNew = 0.0; double sumImaginaryNew = 0.0; @@ -622,6 +744,7 @@ reduction(+:energyRecipNew) } #endif energyRecipOld = sysPotRef.boxEnergy[box].recip; + GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_SWAP_ENERGY); } return energyRecipNew - energyRecipOld; } @@ -629,11 +752,11 @@ reduction(+:energyRecipNew) //calculate reciprocal term for inserting some molecules (kindA) in destination // box and removing molecule (kindB) from destination box -double Ewald::SwapRecip(const std::vector &newMol, - const std::vector &oldMol, - const std::vector molIndexNew, - const std::vector molIndexOld, - bool first_call) +double Ewald::MolExchangeReciprocal(const std::vector &newMol, + const std::vector &oldMol, + const std::vector &molIndexNew, + const std::vector &molIndexOld, + bool first_call) { double energyRecipNew = 0.0; double energyRecipOld = 0.0; @@ -641,6 +764,7 @@ double Ewald::SwapRecip(const std::vector &newMol, uint box = newMol[0].GetBox(); if (box < BOXES_WITH_U_NB) { + GOMC_EVENT_START(1, GomcProfileEvent::RECIP_MEMC_ENERGY); uint lengthNew, lengthOld; MoleculeKind const& thisKindNew = newMol[0].GetKind(); MoleculeKind const& thisKindOld = oldMol[0].GetKind(); @@ -648,17 +772,11 @@ double Ewald::SwapRecip(const std::vector &newMol, lengthOld = thisKindOld.NumAtoms(); #ifdef _OPENMP -#if GCC_VERSION >= 90000 #pragma omp parallel for default(none) shared(box, first_call, lengthNew, lengthOld, \ newMol, oldMol, thisKindNew, thisKindOld, molIndexNew, molIndexOld) \ reduction(+:energyRecipNew) -#else - #pragma omp parallel for default(none) shared(box, first_call, lengthNew, lengthOld, \ - newMol, oldMol, thisKindNew, thisKindOld) \ -reduction(+:energyRecipNew) #endif -#endif - for (int i = 0; i < imageSizeRef[box]; i++) { + for (int i = 0; i < (int) imageSizeRef[box]; i++) { double sumRealNew = 0.0; double sumImaginaryNew = 0.0; @@ -729,8 +847,19 @@ reduction(+:energyRecipNew) // Keep hold of the old recip value energyRecipOld = sysPotRef.boxEnergy[box].recip; + GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_MEMC_ENERGY); } +//Because MolExchangeReciprocal does not have a matching GPU function, this is +//a stub function to copy the CPU sumRnew and sumInew vectors to the GPU in case +//the move is accepted. If this function is ported to the GPU, this call should +//be moved to the beginning of the MolExchangeReciprocal function and called +//instead of running the CPU code. +#ifdef GOMC_CUDA + CallMolExchangeReciprocalGPU(ff.particles->getCUDAVars(), imageSizeRef[box], + sumRnew[box], sumInew[box], box); +#endif + // Return the difference between old and new reciprocal energies return energyRecipNew - energyRecipOld; } @@ -772,12 +901,12 @@ void Ewald::RecipInitOrth(uint box, BoxDimensions const& boxAxes) double alpsqr4 = 1.0 / (4.0 * ff.alphaSq[box]); XYZ constValue = boxAxes.axis.Get(box); constValue.Inverse(); - constValue *= 2 * M_PI; + constValue *= 2.0 * M_PI; - double vol = boxAxes.volume[box] / (4 * M_PI); - nkx_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).x / (2 * M_PI)) + 1; - nky_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).y / (2 * M_PI)) + 1; - nkz_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).z / (2 * M_PI)) + 1; + double vol = boxAxes.volume[box] / (4.0 * M_PI); + nkx_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).x / (2.0 * M_PI)) + 1; + nky_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).y / (2.0 * M_PI)) + 1; + nkz_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).z / (2.0 * M_PI)) + 1; kmax[box] = std::max(std::max(nkx_max, nky_max), std::max(nky_max, nkz_max)); for(x = 0; x <= nkx_max; x++) { @@ -832,12 +961,12 @@ void Ewald::RecipInitNonOrth(uint box, BoxDimensions const& boxAxes) cellB.Scale(2, boxAxes.axis.Get(box).z); XYZArray cellB_Inv(3); double det = cellB.AdjointMatrix(cellB_Inv); - cellB_Inv.ScaleRange(0, 3, (2 * M_PI) / det); + cellB_Inv.ScaleRange(0, 3, (2.0 * M_PI) / det); - double vol = boxAxes.volume[box] / (4 * M_PI); - nkx_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).x / (2 * M_PI)) + 1; - nky_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).y / (2 * M_PI)) + 1; - nkz_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).z / (2 * M_PI)) + 1; + double vol = boxAxes.volume[box] / (4.0 * M_PI); + nkx_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).x / (2.0 * M_PI)) + 1; + nky_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).y / (2.0 * M_PI)) + 1; + nkz_max = int(ff.recip_rcut[box] * boxAxes.axis.Get(box).z / (2.0 * M_PI)) + 1; kmax[box] = std::max(std::max(nkx_max, nky_max), std::max(nky_max, nkz_max)); for (x = 0; x <= nkx_max; x++) { @@ -901,11 +1030,11 @@ void Ewald::RecipCountInit(uint box, BoxDimensions const& boxAxes) cellB.Scale(2, constValue.z); XYZArray cellB_Inv(3); double det = cellB.AdjointMatrix(cellB_Inv); - cellB_Inv.ScaleRange(0, 3, (2 * M_PI) / det); + cellB_Inv.ScaleRange(0, 3, (2.0 * M_PI) / det); - nkx_max = int(ff.recip_rcut[box] * constValue.x / (2 * M_PI)) + 1; - nky_max = int(ff.recip_rcut[box] * constValue.y / (2 * M_PI)) + 1; - nkz_max = int(ff.recip_rcut[box] * constValue.z / (2 * M_PI)) + 1; + nkx_max = int(ff.recip_rcut[box] * constValue.x / (2.0 * M_PI)) + 1; + nky_max = int(ff.recip_rcut[box] * constValue.y / (2.0 * M_PI)) + 1; + nkz_max = int(ff.recip_rcut[box] * constValue.z / (2.0 * M_PI)) + 1; for(x = 0; x <= nkx_max; x++) { if(x == 0.0) @@ -950,12 +1079,9 @@ void Ewald::SetRecipRef(uint box) std::memcpy(prefactRef[box], prefact[box], sizeof(double) *imageSize[box]); } #ifdef GOMC_CUDA - CopyCurrentToRefCUDA(ff.particles->getCUDAVars(), - box, imageSize[box]); + CopyCurrentToRefCUDA(ff.particles->getCUDAVars(), box, imageSize[box]); #endif - for(uint b = 0; b < BOXES_WITH_U_NB; b++) { - imageSizeRef[b] = imageSize[b]; - } + imageSizeRef[box] = imageSize[box]; } //calculate correction term for a molecule, with system lambda @@ -964,6 +1090,7 @@ double Ewald::MolCorrection(uint molIndex, uint box) const if (box >= BOXES_WITH_U_NB) return 0.0; + GOMC_EVENT_START(1, GomcProfileEvent::CORR_MOL); double dist, distSq; double correction = 0.0; XYZ virComponents; @@ -978,14 +1105,21 @@ double Ewald::MolCorrection(uint molIndex, uint box) const continue; } for (uint j = i + 1; j < atomSize; j++) { - currentAxes.InRcut(distSq, virComponents, currentCoords, + if(currentAxes.InRcut(distSq, virComponents, currentCoords, start + i, start + j, box); dist = sqrt(distSq); correction += (thisKind.AtomCharge(i) * thisKind.AtomCharge(j) * erf(ff.alpha[box] * dist) / dist); } + // within cutoff + // CalcCoul = qiqj/r(erfc(alpha*r)+erf(alpha*r)) + + // outside cutoff + // erf(alpha*r)*qiqj/r + // CalcCoul = erf } + GOMC_EVENT_STOP(1, GomcProfileEvent::CORR_MOL); return -1.0 * num::qqFact * correction * lambdaCoef * lambdaCoef; } @@ -1028,11 +1162,12 @@ void Ewald::ChangeCorrection(Energy *energyDiff, Energy &dUdL_Coul, } //calculate self term for a box, using system lambda -double Ewald::BoxSelf(BoxDimensions const& boxAxes, uint box) const +double Ewald::BoxSelf(uint box) const { if (box >= BOXES_WITH_U_NB) return 0.0; + GOMC_EVENT_START(1, GomcProfileEvent::SELF_BOX); double self = 0.0; double molSelfEnergy; uint i, j, length, molNum; @@ -1061,8 +1196,10 @@ double Ewald::BoxSelf(BoxDimensions const& boxAxes, uint box) const } } - self = -1.0 * self * ff.alpha[box] * num::qqFact / sqrt(M_PI); + // M_2_SQRTPI is 2/sqrt(PI), so need to multiply by 0.5 to get sqrt(PI) + self *= -1.0 * ff.alpha[box] * num::qqFact * M_2_SQRTPI * 0.5; + GOMC_EVENT_STOP(1, GomcProfileEvent::SELF_BOX); return self; } @@ -1075,12 +1212,12 @@ Virial Ewald::VirialReciprocal(Virial& virial, uint box) const if (box >= BOXES_WITH_U_NB) return tempVir; + GOMC_EVENT_START(1, GomcProfileEvent::RECIP_BOX_VIRIAL); double wT11 = 0.0, wT12 = 0.0, wT13 = 0.0; double wT22 = 0.0, wT23 = 0.0, wT33 = 0.0; - double recipIntra = 0.0; double constVal = 1.0 / (4.0 * ff.alphaSq[box]); - double charge, lambdaCoef; + double lambdaCoef; uint p, length, startAtom, atom; MoleculeLookup::box_iterator thisMol = molLookup.BoxBegin(box), @@ -1091,7 +1228,7 @@ Virial Ewald::VirialReciprocal(Virial& virial, uint box) const #ifdef GOMC_CUDA int numberOfAtoms = 0, atomIndex = 0; - for(int k = 0; k < mols.GetKindsCount(); k++) { + for(int k = 0; k < (int) mols.GetKindsCount(); k++) { MoleculeKind const& thisKind = mols.kinds[k]; numberOfAtoms += thisKind.NumAtoms() * molLookup.NumKindInBox(k, box); } @@ -1131,7 +1268,7 @@ Virial Ewald::VirialReciprocal(Virial& virial, uint box) const #ifdef _OPENMP #pragma omp parallel for default(none) shared(box, constVal) reduction(+:wT11, wT22, wT33) #endif - for (int i = 0; i < imageSizeRef[box]; i++) { + for (int i = 0; i < (int) imageSizeRef[box]; i++) { double factor = prefactRef[box][i] * (sumRref[box][i] * sumRref[box][i] + sumIref[box][i] * sumIref[box][i]); @@ -1163,12 +1300,12 @@ Virial Ewald::VirialReciprocal(Virial& virial, uint box) const currentAxes.UnwrapPBC(atomC, box, comC); diffC = atomC - comC; //scale the charge with lambda for Free energy calc - charge = particleCharge[atom] * lambdaCoef; + double charge = particleCharge[atom] * lambdaCoef; #ifdef _OPENMP #pragma omp parallel for default(none) shared(atom, box, charge, diffC) reduction(+:wT11, wT22, wT33) #endif - for (int i = 0; i < imageSizeRef[box]; i++) { + for (int i = 0; i < (int) imageSizeRef[box]; i++) { //compute the dot product of k and r double arg = Dot(atom, kxRef[box][i], kyRef[box][i], kzRef[box][i], currentCoords); @@ -1202,6 +1339,7 @@ Virial Ewald::VirialReciprocal(Virial& virial, uint box) const // setting virial of reciprocal space tempVir.recip = wT11 + wT22 + wT33; + GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_BOX_VIRIAL); return tempVir; } @@ -1216,6 +1354,7 @@ double Ewald::SwapCorrection(const cbmc::TrialMol& trialMol) const if (box >= BOXES_WITH_U_NB) return 0.0; + GOMC_EVENT_START(1, GomcProfileEvent::CORR_SWAP); double dist, distSq; double correction = 0.0; XYZ virComponents; @@ -1232,6 +1371,7 @@ double Ewald::SwapCorrection(const cbmc::TrialMol& trialMol) const erf(ff.alpha[box] * dist) / dist); } } + GOMC_EVENT_STOP(1, GomcProfileEvent::CORR_SWAP); return num::qqFact * correction; } @@ -1244,6 +1384,7 @@ double Ewald::SwapCorrection(const cbmc::TrialMol& trialMol, if (box >= BOXES_WITH_U_NB) return 0.0; + GOMC_EVENT_START(1, GomcProfileEvent::CORR_SWAP); double dist, distSq; double correction = 0.0; XYZ virComponents; @@ -1265,18 +1406,20 @@ double Ewald::SwapCorrection(const cbmc::TrialMol& trialMol, erf(ff.alpha[box] * dist) / dist); } } + GOMC_EVENT_STOP(1, GomcProfileEvent::CORR_SWAP); return num::qqFact * correction * lambdaCoef * lambdaCoef; } //It's called if we transfer one molecule from one box to another //No need to scale the charge with lambda, since this function is not being -// called from free energy or CFCMC +// called from free energy or NeMTMC double Ewald::SwapSelf(const cbmc::TrialMol& trialMol) const { uint box = trialMol.GetBox(); if (box >= BOXES_WITH_U_NB) return 0.0; + GOMC_EVENT_START(1, GomcProfileEvent::SELF_SWAP); MoleculeKind const& thisKind = trialMol.GetKind(); uint atomSize = thisKind.NumAtoms(); double en_self = 0.0; @@ -1284,7 +1427,9 @@ double Ewald::SwapSelf(const cbmc::TrialMol& trialMol) const for (uint i = 0; i < atomSize; i++) { en_self -= (thisKind.AtomCharge(i) * thisKind.AtomCharge(i)); } - return (en_self * ff.alpha[box] * num::qqFact / sqrt(M_PI)); + GOMC_EVENT_STOP(1, GomcProfileEvent::SELF_SWAP); + // M_2_SQRTPI is 2/sqrt(PI), so need to multiply by 0.5 to get sqrt(PI) + return (en_self * ff.alpha[box] * num::qqFact * M_2_SQRTPI * 0.5); } //It's called in free energy calculation to calculate the change in @@ -1302,7 +1447,8 @@ void Ewald::ChangeSelf(Energy *energyDiff, Energy &dUdL_Coul, for (uint i = 0; i < atomSize; i++) { en_self += (particleCharge[i + start] * particleCharge[i + start]); } - en_self *= -1.0 * ff.alpha[box] * num::qqFact / sqrt(M_PI); + // M_2_SQRTPI is 2/sqrt(PI), so need to multiply by 0.5 to get sqrt(PI) + en_self *= -1.0 * ff.alpha[box] * num::qqFact * M_2_SQRTPI * 0.5; //Calculate the energy difference for each lambda state for (uint s = 0; s < lambdaSize; s++) { @@ -1331,6 +1477,28 @@ void Ewald::UpdateRecip(uint box) #endif } +//copy reciprocal values from ref to new +//Used to reinitialize these two arrays for MultiParticle moves so that the +//forces can be computed using the current positions, before making the MP +//move. The other option would be to create different versions of these +//functions to access SumIref and SumRref instead of SumInew and SumRnew. +void Ewald::CopyRecip(uint box) +{ + if (box >= BOXES_WITH_U_NB) + return; + +#ifdef _OPENMP + #pragma omp parallel default(none) shared(box) +#endif + { + std::memcpy(sumRnew[box], sumRref[box], sizeof(double) * imageSizeRef[box]); + std::memcpy(sumInew[box], sumIref[box], sizeof(double) * imageSizeRef[box]); + } +#ifdef GOMC_CUDA + CopyRefToNewCUDA(ff.particles->getCUDAVars(), box, imageSizeRef[box]); +#endif +} + void Ewald::UpdateRecipVec(uint box) { double *tempKx, *tempKy, *tempKz, *tempHsqr, *tempPrefact; @@ -1351,13 +1519,12 @@ void Ewald::UpdateRecipVec(uint box) kz[box] = tempKz; hsqr[box] = tempHsqr; prefact[box] = tempPrefact; + + imageSizeRef[box] = imageSize[box]; + #ifdef GOMC_CUDA UpdateRecipVecCUDA(ff.particles->getCUDAVars(), box); #endif - - for(uint b = 0; b < BOXES_WITH_U_NB; b++) { - imageSizeRef[b] = imageSize[b]; - } } void compareDouble(const double &x, const double &y, const int &i) @@ -1375,23 +1542,35 @@ void Ewald::BoxForceReciprocal(XYZArray const& molCoords, uint box) { if(multiParticleEnabled && (box < BOXES_WITH_U_NB)) { - double constValue = 2.0 * ff.alpha[box] / sqrt(M_PI); + GOMC_EVENT_START(1, GomcProfileEvent::RECIP_BOX_FORCE); + // M_2_SQRTPI is 2/sqrt(PI) + double constValue = ff.alpha[box] * M_2_SQRTPI; #ifdef GOMC_CUDA + bool *particleUsed; + particleUsed = new bool[atomForceRec.Count()]; + memset((void *) particleUsed, false, atomForceRec.Count() * sizeof(bool)); +#if ENSEMBLE == GEMC || ENSEMBLE == GCMC + memset((void *) particleUsed, false, atomForceRec.Count() * sizeof(bool)); MoleculeLookup::box_iterator thisMol = molLookup.BoxBegin(box); MoleculeLookup::box_iterator end = molLookup.BoxEnd(box); while(thisMol != end) { uint molIndex = *thisMol; uint length = mols.GetKind(molIndex).NumAtoms(); uint start = mols.MolStart(molIndex); - for(int p = 0; p < length; p++) { + for(uint p = 0; p < length; p++) { atomForceRec.Set(start + p, 0.0, 0.0, 0.0); + particleUsed[start + p] = true; } molForceRec.Set(molIndex, 0.0, 0.0, 0.0); thisMol++; } - // initialize the start and end of each box - initializeBoxRange(); +#else + //Only one box, so clear all atoms and molecules and mark all particles as Used + atomForceRec.Reset(); + molForceRec.Reset(); + memset((void *) particleUsed, true, atomForceRec.Count() * sizeof(bool)); +#endif CallBoxForceReciprocalGPU( ff.particles->getCUDAVars(), @@ -1399,21 +1578,20 @@ void Ewald::BoxForceReciprocal(XYZArray const& molCoords, molForceRec, particleCharge, particleMol, - particleKind, particleHasNoCharge, + particleUsed, startMol, lengthMol, ff.alpha[box], ff.alphaSq[box], num::qqFact, constValue, - imageSize[box], + imageSizeRef[box], molCoords, - boxStart[box], - boxEnd[box], currentAxes, box ); + delete[] particleUsed; #else // molecule iterator MoleculeLookup::box_iterator thisMol = molLookup.BoxBegin(box); @@ -1431,7 +1609,6 @@ void Ewald::BoxForceReciprocal(XYZArray const& molCoords, for(p = start; p < start + length; p++) { double X = 0.0, Y = 0.0, Z = 0.0; - double intraForce; if(!particleHasNoCharge[p]) { // subtract the intra forces(correction) @@ -1442,7 +1619,7 @@ void Ewald::BoxForceReciprocal(XYZArray const& molCoords, double dist = sqrt(distSq); double expConstValue = exp(-1.0 * ff.alphaSq[box] * distSq); double qiqj = particleCharge[p] * particleCharge[j] * num::qqFact; - intraForce = qiqj * lambdaCoef * lambdaCoef / distSq; + double intraForce = qiqj * lambdaCoef * lambdaCoef / distSq; intraForce *= ((erf(ff.alpha[box] * dist) / dist) - constValue * expConstValue); X -= intraForce * distVect.x; @@ -1453,15 +1630,15 @@ void Ewald::BoxForceReciprocal(XYZArray const& molCoords, #ifdef _OPENMP #pragma omp parallel for default(none) shared(box, lambdaCoef, molCoords, p) reduction(+:X, Y, Z) #endif - for(int i = 0; i < imageSize[box]; i++) { - double dot = Dot(p, kx[box][i], ky[box][i], kz[box][i], molCoords); + for(int i = 0; i < (int) imageSizeRef[box]; i++) { + double dot = Dot(p, kxRef[box][i], kyRef[box][i], kzRef[box][i], molCoords); - double factor = 2.0 * particleCharge[p] * prefact[box][i] * lambdaCoef * + double factor = 2.0 * particleCharge[p] * prefactRef[box][i] * lambdaCoef * (sin(dot) * sumRnew[box][i] - cos(dot) * sumInew[box][i]); - X += factor * kx[box][i]; - Y += factor * ky[box][i]; - Z += factor * kz[box][i]; + X += factor * kxRef[box][i]; + Y += factor * kyRef[box][i]; + Z += factor * kzRef[box][i]; } } atomForceRec.Set(p, X, Y, Z); @@ -1470,28 +1647,13 @@ void Ewald::BoxForceReciprocal(XYZArray const& molCoords, thisMol++; } #endif + GOMC_EVENT_STOP(1, GomcProfileEvent::RECIP_BOX_FORCE); } } double Ewald::GetLambdaCoef(uint molA, uint box) const { - double lambda = lambdaRef.GetLambdaCoulomb(molA, mols.GetMolKind(molA), box); + double lambda = lambdaRef.GetLambdaCoulomb(molA, box); //Each charge gets sq root of it. return sqrt(lambda); } - -void Ewald::initializeBoxRange() -{ - if(BOX_TOTAL == 1) { - boxStart[mv::BOX0] = 0; - boxEnd[mv::BOX0] = currentCoords.Count(); - } else { - boxStart[mv::BOX0] = 0; - MoleculeLookup::box_iterator startBox1 = molLookup.BoxBegin(mv::BOX1); - int firstMoleculeInBox1 = *startBox1; - int indexOfFirstAtomInBox1 = mols.MolStart(firstMoleculeInBox1); - boxEnd[mv::BOX0] = indexOfFirstAtomInBox1; - boxStart[mv::BOX1] = indexOfFirstAtomInBox1; - boxEnd[mv::BOX1] = currentCoords.Count(); - } -} From e62e27f6af020fa375f56100edb4e10054836d46 Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Mon, 22 Nov 2021 19:05:08 -0600 Subject: [PATCH 07/16] Revert "Merge branch 'development' of https://github.com/GregorySchwing/GOMC into development" This reverts commit 47383bc5e66ddd447cddce88c141961dc50cdb19, reversing changes made to 6557396e04166d3d41a88a0b806d1035add2a33d. --- src/Ewald.cpp | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/src/Ewald.cpp b/src/Ewald.cpp index 8d19fa151..fa091d8b6 100644 --- a/src/Ewald.cpp +++ b/src/Ewald.cpp @@ -1105,18 +1105,12 @@ double Ewald::MolCorrection(uint molIndex, uint box) const continue; } for (uint j = i + 1; j < atomSize; j++) { - if(currentAxes.InRcut(distSq, virComponents, currentCoords, + currentAxes.InRcut(distSq, virComponents, currentCoords, start + i, start + j, box); dist = sqrt(distSq); correction += (thisKind.AtomCharge(i) * thisKind.AtomCharge(j) * erf(ff.alpha[box] * dist) / dist); } - // within cutoff - // CalcCoul = qiqj/r(erfc(alpha*r)+erf(alpha*r)) - - // outside cutoff - // erf(alpha*r)*qiqj/r - // CalcCoul = erf } GOMC_EVENT_STOP(1, GomcProfileEvent::CORR_MOL); From f5e3ed7878de05a2698913de2fd25ee581e77572 Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Mon, 6 Dec 2021 18:21:13 -0600 Subject: [PATCH 08/16] add tests to dev --- test/Run_Examples.py | 195 ++++++++++++++++++++++++++++++----------- test/Setup_Examples.sh | 28 +++--- test/runIntTests.sh | 24 +++++ 3 files changed, 181 insertions(+), 66 deletions(-) create mode 100755 test/runIntTests.sh diff --git a/test/Run_Examples.py b/test/Run_Examples.py index fc4aa1b05..925f751d9 100644 --- a/test/Run_Examples.py +++ b/test/Run_Examples.py @@ -1,5 +1,14 @@ #!/usr/bin/python +HEADER = '\033[95m' +OKBLUE = '\033[94m' +OKGREEN = '\033[92m' +WARNING = '\033[93m' +FAIL = '\033[91m' +ENDC = '\033[0m' +BOLD = '\033[1m' +UNDERLINE = '\033[4m' + def substring_after(s, delim): return s.partition(delim)[2] @@ -11,10 +20,11 @@ def substring_after(s, delim): import datetime import subprocess from subprocess import Popen, PIPE, STDOUT +from filecmp import cmp binaries_dict = {} GPU_binaries_dict = {} -os.chdir("new_feature_binaries") +os.chdir("new_binaries") pathToDir = os.getcwd() binaries_cpu_new = sorted(glob.glob('GOMC_CPU_*'), key=os.path.getmtime) binaries_gpu_new = sorted(glob.glob('GOMC_GPU_*'), key=os.path.getmtime) @@ -52,72 +62,155 @@ def substring_after(s, delim): # print(len(path) * '---', file) if file==confFileName: newOrRef = "" - if "new_feature_cpu" in path: + cpuOrGpu = "" + if "new_cpu" in path: newOrRef = "_new" + cpuOrGpu = "CPU_" elif "ref_cpu" in path: newOrRef = "_ref" + cpuOrGpu = "CPU_" + elif "new_gpu" in path: + newOrRef = "_new" + cpuOrGpu = "GPU_" + + elif "ref_gpu" in path: + newOrRef = "_ref" + cpuOrGpu = "GPU_" - run = False - if cpuOrGpu+"NVT"+newOrRef in binaries_dict and 'NVT' in path and 'GEMC_NVT' not in path: + if cpuOrGpu+"NVT"+newOrRef in binaries_dict and 'NVT' in path and 'NVT_GEMC' not in path: print("Call GOMC") - command = binaries_dict[cpuOrGpu+"NVT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NVT"+newOrRef,os.path.basename(root) - print(binaries_dict[cpuOrGpu+"NVT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NVT"+newOrRef,os.path.basename(root)) + command = binaries_dict[cpuOrGpu+"NVT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NVT"+newOrRef,cpuOrGpu,newOrRef,"NVT_"+os.path.basename(root) + print(binaries_dict[cpuOrGpu+"NVT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NVT"+newOrRef,cpuOrGpu,newOrRef,"NVT_"+os.path.basename(root)) listOfTests.append(command) - run = True - elif cpuOrGpu+"NPT"+newOrRef in binaries_dict and 'NPT' in path and 'GEMC_NPT' not in path: + elif cpuOrGpu+"NPT"+newOrRef in binaries_dict and 'NPT' in path and 'NPT_GEMC' not in path: print("Call GOMC") - print(binaries_dict[cpuOrGpu+"NPT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NPT"+newOrRef,os.path.basename(root)) - command = binaries_dict[cpuOrGpu+"NPT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NPT"+newOrRef,os.path.basename(root) + print(binaries_dict[cpuOrGpu+"NPT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NPT"+newOrRef,cpuOrGpu,newOrRef,"NPT_"+os.path.basename(root)) + command = binaries_dict[cpuOrGpu+"NPT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NPT"+newOrRef,cpuOrGpu,newOrRef,"NPT_"+os.path.basename(root) listOfTests.append(command) - run = True elif cpuOrGpu+"GCMC"+newOrRef in binaries_dict and 'GCMC' in path: print("Call GOMC") - print(binaries_dict[cpuOrGpu+"GCMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GCMC"+newOrRef,os.path.basename(root)) - command = binaries_dict[cpuOrGpu+"GCMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GCMC"+newOrRef,os.path.basename(root) + print(binaries_dict[cpuOrGpu+"GCMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GCMC"+newOrRef,cpuOrGpu,newOrRef,"GCMC_"+os.path.basename(root)) + command = binaries_dict[cpuOrGpu+"GCMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GCMC"+newOrRef,cpuOrGpu,newOrRef,"GCMC_"+os.path.basename(root) listOfTests.append(command) - run = True - elif cpuOrGpu+"GEMC"+newOrRef in binaries_dict and 'GEMC_NVT' in path: + elif cpuOrGpu+"GEMC"+newOrRef in binaries_dict and 'NVT_GEMC' in path: print("Call GOMC") - print(binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GEMC_NVT"+newOrRef,os.path.basename(root)) - command = binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GEMC_NVT"+newOrRef,os.path.basename(root) + print(binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"NVT_GEMC"+newOrRef,cpuOrGpu,newOrRef,"NPT_GEMC_"+os.path.basename(root)) + command = binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"NVT_GEMC"+newOrRef,cpuOrGpu,newOrRef,"NVT_GEMC_"+os.path.basename(root) listOfTests.append(command) - run = True - elif cpuOrGpu+"GEMC"+newOrRef in binaries_dict and 'GEMC_NPT' in path: + elif cpuOrGpu+"GEMC"+newOrRef in binaries_dict and 'NPT_GEMC' in path: print("Call GOMC") - print(binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GEMC_NPT"+newOrRef,os.path.basename(root)) - command = binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"GEMC_NPT"+newOrRef,os.path.basename(root) + print(binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"NPT_GEMC"+newOrRef,cpuOrGpu,newOrRef,"NPT_GEMC_"+os.path.basename(root)) + command = binaries_dict[cpuOrGpu+"GEMC"+newOrRef],os.path.abspath(root),cpuOrGpu+"NPT_GEMC"+newOrRef,cpuOrGpu,newOrRef,"NPT_GEMC_"+os.path.basename(root) listOfTests.append(command) - run = True - else: - run = False - print(listOfTests) # Create the pandas DataFrame -df = pd.DataFrame(listOfTests, columns = ['PathToBinary', 'PathToExample', 'Binary', 'Example']) +colNames = ['PathToBinary', 'PathToExample', 'Binary', 'CPU_or_GPU','New_or_Ref','Example'] +df = pd.DataFrame(listOfTests, columns = colNames) df = df.sort_values(by=['Example']) print(df) -for index, row in df.iterrows(): - print(index) - print(row) -""" - os.chdir(row['PathToExample']) - write_log_data = "Changing directory to {}\n".format(row['PathToExample']) - Log_Template_file.write(str(write_log_data)) - command = (row['PathToBinary'] + " in.conf > out.log") - write_log_data = "Issuing command: {}\n".format(command) - Log_Template_file.write(str(write_log_data)) - start = time.time() - exec_GOMC_run_command = subprocess.Popen(command, shell=True, stderr=STDOUT) - write_log_data = "Waiting for GOMC Example {} {} to finish.\n".format(row['Binary'],row['Example']) - print(str(write_log_data)) - Log_Template_file.write(str(write_log_data)) - GOMC_pid_status = os.waitpid(exec_GOMC_run_command.pid, os.WSTOPPED) # pauses python until box 0 sim done - end = time.time() - write_log_data = "Elapsed time: {}.\n".format(datetime.timedelta(seconds=end-start)) - Log_Template_file.write(str(write_log_data)) - print(str(write_log_data)) - write_log_data = "The GOMC Example {} {} has finished.\n".format(row['Binary'],row['Example']) - print(str(write_log_data)) - Log_Template_file.write(str(write_log_data)) - Log_Template_file.flush() -""" +#for index, row in df.iterrows(): +# print(index) +# print(row) +grouped = df.groupby(['Example']) +all_examples = df['Example'].unique() +CPU_v_GPU_global = True +New_v_Ref_global = True +cross_bool_global = True +CPU_v_GPU_exists_global = False +New_v_Ref_exists_global = False +cross_exists_global = False +for example in all_examples: + ex_df = grouped.get_group(example) + for index, row in ex_df.iterrows(): + #print("run file {}".format(row['PathToExample']+" conf")) + os.chdir(row['PathToExample']) + write_log_data = "Changing directory to {}\n".format(row['PathToExample']) + Log_Template_file.write(str(write_log_data)) + command = (row['PathToBinary'] + " in.conf > out.log") + write_log_data = "Issuing command: {}\n".format(command) + Log_Template_file.write(str(write_log_data)) + start = time.time() + exec_GOMC_run_command = subprocess.Popen(command, shell=True, stderr=STDOUT) + write_log_data = "Waiting for GOMC Example {} {} to finish.\n".format(row['Binary'],row['Example']) + print(str(write_log_data)) + Log_Template_file.write(str(write_log_data)) + GOMC_pid_status = os.waitpid(exec_GOMC_run_command.pid, os.WSTOPPED) # pauses python until box 0 sim done + end = time.time() + write_log_data = "Elapsed time: {}.\n".format(datetime.timedelta(seconds=end-start)) + Log_Template_file.write(str(write_log_data)) + print(str(write_log_data)) + write_log_data = "The GOMC Example {} {} has finished.\n".format(row['Binary'],row['Example']) + print(str(write_log_data)) + Log_Template_file.write(str(write_log_data)) + Log_Template_file.flush() + + # Create a list of the PDB files in this example + full_path_pdb_files = sorted(glob.glob(os.path.join(row['PathToExample'],'*.pdb')), key=os.path.getmtime) + just_file_names = [] + for path in full_path_pdb_files: + just_file_names.append(os.path.basename(path)) + print(just_file_names) + cross = ex_df.merge(ex_df, on=['Example'],how='outer') + print('---', ex_df['Example'].iloc[0]) + CPU_v_GPU = True + New_v_Ref = True + cross_bool = True + CPU_v_GPU_exists = False + New_v_Ref_exists = False + cross_exists = False + for pdb_file in just_file_names: + print(2 * '---', pdb_file) + my_tuples = [] + for index, row in cross.iterrows(): + f1 = os.path.join(row['PathToExample_x'],pdb_file) + f2 = os.path.join(row['PathToExample_y'],pdb_file) + if ((row['CPU_or_GPU_x'] != row['CPU_or_GPU_y']) and (row['New_or_Ref_x'] == row['New_or_Ref_y'])): + CPU_v_GPU_exists = True + CPU_v_GPU_exists_global = True + result = cmp(f1, f2, shallow=False) + CPU_v_GPU = CPU_v_GPU and result + CPU_v_GPU_global = CPU_v_GPU_global and result + elif ((row['CPU_or_GPU_x'] == row['CPU_or_GPU_y']) and (row['New_or_Ref_x'] != row['New_or_Ref_y'])): + New_v_Ref_exists = True + New_v_Ref_exists_global = True + result = cmp(f1, f2, shallow=False) + New_v_Ref = New_v_Ref and result + New_v_Ref_global = New_v_Ref_global and result + elif ((row['CPU_or_GPU_x'] != row['CPU_or_GPU_y']) and (row['New_or_Ref_x'] != row['New_or_Ref_y'])): + cross_exists = True + cross_exists_global = True + result = cmp(f1, f2, shallow=False) + cross_bool = cross_bool and result + cross_bool_global = cross_bool_global and result + if(CPU_v_GPU_exists): + if(CPU_v_GPU): + print((3 * '---')+"CPU_v_GPU: "+ OKGREEN + "PASS" + ENDC) + else: + print((3 * '---')+"CPU_v_GPU: "+ FAIL + "FAIL" + ENDC) + if(New_v_Ref_exists): + if(New_v_Ref): + print((3 * '---')+"New vs Ref: "+ OKGREEN + "PASS" + ENDC) + else: + print((3 * '---')+"New vs Ref: "+ FAIL + "FAIL" + ENDC) + if(cross_exists): + if(cross_bool): + print((3 * '---')+"CPU vs GPU X New vs Ref: "+ OKGREEN + "PASS" + ENDC) + else: + print((3 * '---')+"CPU vs GPU X New vs Ref: "+ FAIL + "FAIL" + ENDC) + + +if(CPU_v_GPU_exists_global): + if(CPU_v_GPU_global): + print("CPU_v_GPU Global: "+ OKGREEN + "PASS" + ENDC) + else: + print("CPU_v_GPU Global: "+ FAIL + "FAIL" + ENDC) +if(New_v_Ref_exists_global): + if(New_v_Ref_global): + print("New vs Ref Global: "+ OKGREEN + "PASS" + ENDC) + else: + print("New vs Ref Global: "+ FAIL + "FAIL" + ENDC) +if(cross_exists_global): + if(cross_bool_global): + print("CPU vs GPU X New vs Ref Global: "+ OKGREEN + "PASS" + ENDC) + else: + print("CPU vs GPU X New vs Ref Global: "+ FAIL + "FAIL" + ENDC) diff --git a/test/Setup_Examples.sh b/test/Setup_Examples.sh index d5859e467..46c73272c 100755 --- a/test/Setup_Examples.sh +++ b/test/Setup_Examples.sh @@ -1,30 +1,28 @@ #!/bin/bash git clone https://github.com/GOMC-WSU/GOMC_Examples.git mkdir integration -cd integration -mkdir new_cpu -cp -frd GOMC_Examples new_cpu +mkdir integration/new_cpu +cp -frd GOMC_Examples integration/new_cpu -mkdir ref_cpu -cp -frd GOMC_Examples ref_cpu +mkdir integration/ref_cpu +cp -frd GOMC_Examples integration/ref_cpu -mkdir new_gpu -cp -frd GOMC_Examples new_gpu +mkdir integration/new_gpu +cp -frd GOMC_Examples integration/new_gpu -mkdir ref_gpu -cp -frd GOMC_Examples ref_gpu +mkdir integration/ref_gpu +cp -frd GOMC_Examples integration/ref_gpu +cd .. -cd ../.. startingBranch=$(git symbolic-ref HEAD | sed -e 's,.*/\(.*\),\1,') echo "Building $startingBranch binaries" -./metamake.sh -g NVT +./metamake.sh -g NVT NPT GCMC GEMC mkdir -p test/new_binaries cp -frdp ./bin/* test/new_binaries - echo "$startingBranch" if [ $startingBranch == "development" ]; then echo "I am on development; checking out main" @@ -36,9 +34,9 @@ fi rm -frd bin refBranch=$(git symbolic-ref HEAD | sed -e 's,.*/\(.*\),\1,') -echo "Building $refBbranch binaries" -./metamake.sh -g NVT +echo "Building $refBranch binaries" +./metamake.sh -g NVT NPT GCMC GEMC mkdir -p test/ref_binaries cp -frd ./bin/* test/ref_binaries cd test - +git checkout $startingBranch diff --git a/test/runIntTests.sh b/test/runIntTests.sh new file mode 100755 index 000000000..64d1b842f --- /dev/null +++ b/test/runIntTests.sh @@ -0,0 +1,24 @@ +#!/bin/bash +# Example with 28 cores for OpenMP +# +# Project/Account +#SBATCH -A greg +# +# Number of cores +#SBATCH -c 6 -w, --nodelist=potoff33 +# +# Runtime of this jobs is less then 12 hours. +#SBATCH --time=168:00:00 +# +#SBATCH --mail-user=go2432@wayne.edu + +#SBATCH -o output_%j.out + +#SBATCH -e errors_%j.err + + +# Clear the environment from any previously loaded modules +cd /home6/greg/GOMC/test +bash ./Setup_Examples.sh +python Run_Examples.py > integrationTest.log +# End of submit file From 77e414fb4c20730a9290b1430a27cb74701460e8 Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Tue, 7 Dec 2021 11:14:37 -0600 Subject: [PATCH 09/16] build all bins --- test/Setup_Examples.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/Setup_Examples.sh b/test/Setup_Examples.sh index 46c73272c..fe6b94120 100755 --- a/test/Setup_Examples.sh +++ b/test/Setup_Examples.sh @@ -19,7 +19,7 @@ cd .. startingBranch=$(git symbolic-ref HEAD | sed -e 's,.*/\(.*\),\1,') echo "Building $startingBranch binaries" -./metamake.sh -g NVT NPT GCMC GEMC +./metamake.sh -g mkdir -p test/new_binaries cp -frdp ./bin/* test/new_binaries @@ -35,7 +35,7 @@ fi rm -frd bin refBranch=$(git symbolic-ref HEAD | sed -e 's,.*/\(.*\),\1,') echo "Building $refBranch binaries" -./metamake.sh -g NVT NPT GCMC GEMC +./metamake.sh -g mkdir -p test/ref_binaries cp -frd ./bin/* test/ref_binaries cd test From 98bdc7e798c2ceb215018fc1530aae106cea3671 Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Thu, 9 Dec 2021 16:27:23 -0600 Subject: [PATCH 10/16] Fix HVAP. This also fixes undefined behavior from using uninitialized variables --- src/OutputVars.cpp | 50 ++++++++++++++++++++++++++++++---------------- 1 file changed, 33 insertions(+), 17 deletions(-) diff --git a/src/OutputVars.cpp b/src/OutputVars.cpp index 0c9010712..ab56f556c 100644 --- a/src/OutputVars.cpp +++ b/src/OutputVars.cpp @@ -18,6 +18,19 @@ OutputVars::OutputVars(System & sys, StaticVals const& statV, const std::vector< T_in_K(statV.forcefield.T_in_K), calc(sys.calcEnergy), molKindNames(molKindNames) { InitRef(sys, statV); + for (int b = 0; b < BOX_TOTAL; ++b){ + compressability[b] = 0.0; + enthalpy[b] = 0.0; + } + #if ENSEMBLE == GEMC + liqBox = 0; + vapBox = 0; + heatOfVap = 0.0; + for (int b = 0; b < BOX_TOTAL; ++b){ + heatOfVap_energy_term_box[b] = 0.0; + heatOfVap_density_term_box[b] = 0.0; + } + #endif } void OutputVars::InitRef(System & sys, StaticVals const& statV) @@ -101,11 +114,28 @@ void OutputVars::CalcAndConvert(ulong step) molLookupRef->TotalAndDensity(numByBox, numByKindBox, molFractionByKindBox, densityByKindBox, volInvRef); + for (uint b = 0; b < BOX_TOTAL; b++) { + densityTot[b] = 0.0; + for (uint k = 0; k < numKinds; k++) { + double density = densityByKindBox[k + numKinds * b]; + + // Convert density to g/ml (which is equivalent to g/cm3) + // To get kg/m3, multiply output densities by 1000. + density *= unit::MOLECULES_PER_A3_TO_MOL_PER_CM3 * + kindsRef[k].molMass; + densityTot[b] += density; + } + densityTot[b] *= 1000; + } + #if ENSEMBLE == GEMC //Determine which box is liquid for purposes of heat of vap. - if (densityByKindBox[numKinds] > densityByKindBox[0]) { + if (densityTot[mv::BOX1] >= densityTot[mv::BOX0]) { vapBox = mv::BOX0; liqBox = mv::BOX1; + } else { + vapBox = mv::BOX1; + liqBox = mv::BOX0; } #endif @@ -163,8 +193,8 @@ void OutputVars::CalcAndConvert(ulong step) compressability[b] = (pressure[b]) * (volumeRef[b]) / numByBox[b] / (T_in_K) / (UNIT_CONST_H::unit::K_MOLECULE_PER_A3_TO_BAR); enthalpy[b] = (energyRef[b].total / numByBox[b] + rawPressure[b] * volumeRef[b] / numByBox[b]) * UNIT_CONST_H::unit::K_TO_KJ_PER_MOL; } else { - compressability[b] = 0; - enthalpy[b] = 0; + compressability[b] = 0.0; + enthalpy[b] = 0.0; } #if ENSEMBLE == GEMC // delta Hv = (Uv-Ul) + P(Vv-Vl) @@ -194,18 +224,4 @@ void OutputVars::CalcAndConvert(ulong step) } } } - - for (uint b = 0; b < BOX_TOTAL; b++) { - densityTot[b] = 0.0; - for (uint k = 0; k < numKinds; k++) { - double density = densityByKindBox[k + numKinds * b]; - - // Convert density to g/ml (which is equivalent to g/cm3) - // To get kg/m3, multiply output densities by 1000. - density *= unit::MOLECULES_PER_A3_TO_MOL_PER_CM3 * - kindsRef[k].molMass; - densityTot[b] += density; - } - densityTot[b] *= 1000; - } } From 8b6e6b5638016db9599ffdfda293143b32d11069 Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Thu, 9 Dec 2021 17:07:50 -0600 Subject: [PATCH 11/16] compressability -> compressibility spell fix --- src/BlockOutput.cpp | 2 +- src/OutputVars.cpp | 6 +++--- src/OutputVars.h | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/BlockOutput.cpp b/src/BlockOutput.cpp index 448111fb6..e592fcb4e 100644 --- a/src/BlockOutput.cpp +++ b/src/BlockOutput.cpp @@ -217,7 +217,7 @@ void BlockAverages::InitWatchSingle(config_setup::TrackedVars const& tracked) blocks[out::PRESSURE_IDX].SetRef(&var->pressure[b], b); blocks[out::MOL_NUM_IDX].SetRef(&var->numByBox[b], b); blocks[out::DENSITY_IDX].SetRef(&var->densityTot[b], b); - blocks[out::COMPRESSIBILITY_IDX].SetRef(&var->compressability[b], b); + blocks[out::COMPRESSIBILITY_IDX].SetRef(&var->compressibility[b], b); blocks[out::ENTHALPY_IDX].SetRef(&var->enthalpy[b], b); blocks[out::SURF_TENSION_IDX].SetRef(&var->surfaceTens[b], b); #if ENSEMBLE == GEMC diff --git a/src/OutputVars.cpp b/src/OutputVars.cpp index ab56f556c..ed78924a0 100644 --- a/src/OutputVars.cpp +++ b/src/OutputVars.cpp @@ -19,7 +19,7 @@ OutputVars::OutputVars(System & sys, StaticVals const& statV, const std::vector< { InitRef(sys, statV); for (int b = 0; b < BOX_TOTAL; ++b){ - compressability[b] = 0.0; + compressibility[b] = 0.0; enthalpy[b] = 0.0; } #if ENSEMBLE == GEMC @@ -190,10 +190,10 @@ void OutputVars::CalcAndConvert(ulong step) pressure[b] = rawPressure[b]; pressure[b] *= unit::K_MOLECULE_PER_A3_TO_BAR; if(numByBox[b] != 0){ - compressability[b] = (pressure[b]) * (volumeRef[b]) / numByBox[b] / (T_in_K) / (UNIT_CONST_H::unit::K_MOLECULE_PER_A3_TO_BAR); + compressibility[b] = (pressure[b]) * (volumeRef[b]) / numByBox[b] / (T_in_K) / (UNIT_CONST_H::unit::K_MOLECULE_PER_A3_TO_BAR); enthalpy[b] = (energyRef[b].total / numByBox[b] + rawPressure[b] * volumeRef[b] / numByBox[b]) * UNIT_CONST_H::unit::K_TO_KJ_PER_MOL; } else { - compressability[b] = 0.0; + compressibility[b] = 0.0; enthalpy[b] = 0.0; } #if ENSEMBLE == GEMC diff --git a/src/OutputVars.h b/src/OutputVars.h index e5f377451..e56e97f5c 100644 --- a/src/OutputVars.h +++ b/src/OutputVars.h @@ -40,7 +40,7 @@ class OutputVars //Intermediate vars. uint * numByBox, * numByKindBox; double *molFractionByKindBox, *densityByKindBox, - pressure[BOXES_WITH_U_NB], densityTot[BOX_TOTAL], compressability[BOXES_WITH_U_NB], enthalpy[BOXES_WITH_U_NB]; + pressure[BOXES_WITH_U_NB], densityTot[BOX_TOTAL], compressibility[BOXES_WITH_U_NB], enthalpy[BOXES_WITH_U_NB]; double pressureTens[BOXES_WITH_U_NB][3][3]; double surfaceTens[BOXES_WITH_U_NB]; ulong pCalcFreq; From d406de51dea3db9c51ad745dbca4fc1286aa2573 Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Fri, 10 Dec 2021 09:19:59 -0600 Subject: [PATCH 12/16] Fix oob memory in initialization --- src/OutputVars.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OutputVars.cpp b/src/OutputVars.cpp index ed78924a0..6ca5a80d0 100644 --- a/src/OutputVars.cpp +++ b/src/OutputVars.cpp @@ -18,7 +18,7 @@ OutputVars::OutputVars(System & sys, StaticVals const& statV, const std::vector< T_in_K(statV.forcefield.T_in_K), calc(sys.calcEnergy), molKindNames(molKindNames) { InitRef(sys, statV); - for (int b = 0; b < BOX_TOTAL; ++b){ + for (int b = 0; b < BOXES_WITH_U_NB; ++b){ compressibility[b] = 0.0; enthalpy[b] = 0.0; } From 8bea158c25961fc0fcc909150a8f08c7a51c82e0 Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Fri, 10 Dec 2021 09:35:01 -0600 Subject: [PATCH 13/16] Codacy for tests --- test/Run_Examples.py | 24 +++++++++++------------- test/Setup_Examples.sh | 4 ++-- 2 files changed, 13 insertions(+), 15 deletions(-) diff --git a/test/Run_Examples.py b/test/Run_Examples.py index 925f751d9..fe5d38065 100644 --- a/test/Run_Examples.py +++ b/test/Run_Examples.py @@ -14,12 +14,11 @@ def substring_after(s, delim): import os import glob -import sys import pandas as pd import time import datetime import subprocess -from subprocess import Popen, PIPE, STDOUT +from subprocess import PIPE, STDOUT from filecmp import cmp binaries_dict = {} GPU_binaries_dict = {} @@ -72,17 +71,16 @@ def substring_after(s, delim): elif "new_gpu" in path: newOrRef = "_new" cpuOrGpu = "GPU_" - elif "ref_gpu" in path: newOrRef = "_ref" cpuOrGpu = "GPU_" - if cpuOrGpu+"NVT"+newOrRef in binaries_dict and 'NVT' in path and 'NVT_GEMC' not in path: + if cpuOrGpu+"NVT"+newOrRef in binaries_dict and 'NVT' in path and 'NVT_GEMC' not in path: print("Call GOMC") command = binaries_dict[cpuOrGpu+"NVT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NVT"+newOrRef,cpuOrGpu,newOrRef,"NVT_"+os.path.basename(root) print(binaries_dict[cpuOrGpu+"NVT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NVT"+newOrRef,cpuOrGpu,newOrRef,"NVT_"+os.path.basename(root)) listOfTests.append(command) - elif cpuOrGpu+"NPT"+newOrRef in binaries_dict and 'NPT' in path and 'NPT_GEMC' not in path: + elif cpuOrGpu+"NPT"+newOrRef in binaries_dict and 'NPT' in path and 'NPT_GEMC' not in path: print("Call GOMC") print(binaries_dict[cpuOrGpu+"NPT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NPT"+newOrRef,cpuOrGpu,newOrRef,"NPT_"+os.path.basename(root)) command = binaries_dict[cpuOrGpu+"NPT"+newOrRef],os.path.abspath(root),cpuOrGpu+"NPT"+newOrRef,cpuOrGpu,newOrRef,"NPT_"+os.path.basename(root) @@ -130,11 +128,11 @@ def substring_after(s, delim): write_log_data = "Issuing command: {}\n".format(command) Log_Template_file.write(str(write_log_data)) start = time.time() - exec_GOMC_run_command = subprocess.Popen(command, shell=True, stderr=STDOUT) + exec_GOMC_run_command = subprocess.Popen(command, stderr=STDOUT) write_log_data = "Waiting for GOMC Example {} {} to finish.\n".format(row['Binary'],row['Example']) print(str(write_log_data)) Log_Template_file.write(str(write_log_data)) - GOMC_pid_status = os.waitpid(exec_GOMC_run_command.pid, os.WSTOPPED) # pauses python until box 0 sim done + GOMC_pid_status = os.waitpid(exec_GOMC_run_command.pid, os.WSTOPPED) # pauses python until box 0 sim done end = time.time() write_log_data = "Elapsed time: {}.\n".format(datetime.timedelta(seconds=end-start)) Log_Template_file.write(str(write_log_data)) @@ -148,8 +146,8 @@ def substring_after(s, delim): full_path_pdb_files = sorted(glob.glob(os.path.join(row['PathToExample'],'*.pdb')), key=os.path.getmtime) just_file_names = [] for path in full_path_pdb_files: - just_file_names.append(os.path.basename(path)) - print(just_file_names) + just_file_names.append(os.path.basename(path)) + print(just_file_names) cross = ex_df.merge(ex_df, on=['Example'],how='outer') print('---', ex_df['Example'].iloc[0]) CPU_v_GPU = True @@ -174,14 +172,14 @@ def substring_after(s, delim): New_v_Ref_exists = True New_v_Ref_exists_global = True result = cmp(f1, f2, shallow=False) - New_v_Ref = New_v_Ref and result - New_v_Ref_global = New_v_Ref_global and result + New_v_Ref = New_v_Ref and result + New_v_Ref_global = New_v_Ref_global and result elif ((row['CPU_or_GPU_x'] != row['CPU_or_GPU_y']) and (row['New_or_Ref_x'] != row['New_or_Ref_y'])): cross_exists = True cross_exists_global = True result = cmp(f1, f2, shallow=False) - cross_bool = cross_bool and result - cross_bool_global = cross_bool_global and result + cross_bool = cross_bool and result + cross_bool_global = cross_bool_global and result if(CPU_v_GPU_exists): if(CPU_v_GPU): print((3 * '---')+"CPU_v_GPU: "+ OKGREEN + "PASS" + ENDC) diff --git a/test/Setup_Examples.sh b/test/Setup_Examples.sh index fe6b94120..53491e349 100755 --- a/test/Setup_Examples.sh +++ b/test/Setup_Examples.sh @@ -24,7 +24,7 @@ mkdir -p test/new_binaries cp -frdp ./bin/* test/new_binaries echo "$startingBranch" -if [ $startingBranch == "development" ]; then +if [ "$startingBranch" == "development" ]; then echo "I am on development; checking out main" git checkout main else @@ -39,4 +39,4 @@ echo "Building $refBranch binaries" mkdir -p test/ref_binaries cp -frd ./bin/* test/ref_binaries cd test -git checkout $startingBranch +git checkout "$startingBranch" From 64c68528543241162c5e958d229dccbc6da01ac1 Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Fri, 10 Dec 2021 09:49:48 -0600 Subject: [PATCH 14/16] Fix test codacy --- test/Run_Examples.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/Run_Examples.py b/test/Run_Examples.py index fe5d38065..cb97efe2e 100644 --- a/test/Run_Examples.py +++ b/test/Run_Examples.py @@ -18,7 +18,7 @@ def substring_after(s, delim): import time import datetime import subprocess -from subprocess import PIPE, STDOUT +from subprocess import STDOUT from filecmp import cmp binaries_dict = {} GPU_binaries_dict = {} @@ -141,7 +141,7 @@ def substring_after(s, delim): print(str(write_log_data)) Log_Template_file.write(str(write_log_data)) Log_Template_file.flush() - + # Create a list of the PDB files in this example full_path_pdb_files = sorted(glob.glob(os.path.join(row['PathToExample'],'*.pdb')), key=os.path.getmtime) just_file_names = [] From e9292ae790bd2fdf047cf95811897e655016d89b Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Fri, 10 Dec 2021 09:57:11 -0600 Subject: [PATCH 15/16] consolidate GEMC pragmas in Output Vars and move outside the loop of boxes so we dont calculate 2x HVAP --- src/OutputVars.cpp | 73 +++++++++++++++++++++++----------------------- 1 file changed, 37 insertions(+), 36 deletions(-) diff --git a/src/OutputVars.cpp b/src/OutputVars.cpp index 6ca5a80d0..44cabdd0e 100644 --- a/src/OutputVars.cpp +++ b/src/OutputVars.cpp @@ -26,7 +26,7 @@ OutputVars::OutputVars(System & sys, StaticVals const& statV, const std::vector< liqBox = 0; vapBox = 0; heatOfVap = 0.0; - for (int b = 0; b < BOX_TOTAL; ++b){ + for (int b = 0; b < BOXES_WITH_U_NB; ++b){ heatOfVap_energy_term_box[b] = 0.0; heatOfVap_density_term_box[b] = 0.0; } @@ -128,17 +128,6 @@ void OutputVars::CalcAndConvert(ulong step) densityTot[b] *= 1000; } -#if ENSEMBLE == GEMC - //Determine which box is liquid for purposes of heat of vap. - if (densityTot[mv::BOX1] >= densityTot[mv::BOX0]) { - vapBox = mv::BOX0; - liqBox = mv::BOX1; - } else { - vapBox = mv::BOX1; - liqBox = mv::BOX0; - } -#endif - for (uint b = 0; b < BOXES_WITH_U_NB; b++) { //Account for dimensionality of virial (raw "virial" is actually a //multiple of the true virial, based on the dimensions stress is exerted @@ -196,32 +185,44 @@ void OutputVars::CalcAndConvert(ulong step) compressibility[b] = 0.0; enthalpy[b] = 0.0; } + } + } + } #if ENSEMBLE == GEMC - // delta Hv = (Uv-Ul) + P(Vv-Vl) - if (numByBox[vapBox] != 0){ - heatOfVap_energy_term_box[vapBox] = energyRef[vapBox].total / numByBox[vapBox]; - heatOfVap_density_term_box[vapBox] = volumeRef[vapBox] / numByBox[vapBox]; - } else { - heatOfVap_energy_term_box[vapBox] = 0.0; - heatOfVap_density_term_box[vapBox] = 0.0; - } - - if (numByBox[liqBox] != 0){ - heatOfVap_energy_term_box[liqBox] = energyRef[liqBox].total / numByBox[liqBox]; - heatOfVap_density_term_box[liqBox] = volumeRef[liqBox] / numByBox[liqBox]; - } else { - heatOfVap_energy_term_box[liqBox] = 0.0; - heatOfVap_density_term_box[liqBox] = 0.0; - } - - heatOfVap = heatOfVap_energy_term_box[vapBox] - - heatOfVap_energy_term_box[liqBox]; - heatOfVap += rawPressure[vapBox] * - (heatOfVap_density_term_box[vapBox] - - heatOfVap_density_term_box[liqBox]); - heatOfVap *= unit::K_TO_KJ_PER_MOL; -#endif + if (pressureCalc) { + if((step + 1) % pCalcFreq == 0 || step == 0) { + //Determine which box is liquid for purposes of heat of vap. + if (densityTot[mv::BOX1] >= densityTot[mv::BOX0]) { + vapBox = mv::BOX0; + liqBox = mv::BOX1; + } else { + vapBox = mv::BOX1; + liqBox = mv::BOX0; + } + // delta Hv = (Uv-Ul) + P(Vv-Vl) + if (numByBox[vapBox] != 0){ + heatOfVap_energy_term_box[vapBox] = energyRef[vapBox].total / numByBox[vapBox]; + heatOfVap_density_term_box[vapBox] = volumeRef[vapBox] / numByBox[vapBox]; + } else { + heatOfVap_energy_term_box[vapBox] = 0.0; + heatOfVap_density_term_box[vapBox] = 0.0; } + + if (numByBox[liqBox] != 0){ + heatOfVap_energy_term_box[liqBox] = energyRef[liqBox].total / numByBox[liqBox]; + heatOfVap_density_term_box[liqBox] = volumeRef[liqBox] / numByBox[liqBox]; + } else { + heatOfVap_energy_term_box[liqBox] = 0.0; + heatOfVap_density_term_box[liqBox] = 0.0; + } + + heatOfVap = heatOfVap_energy_term_box[vapBox] - + heatOfVap_energy_term_box[liqBox]; + heatOfVap += rawPressure[vapBox] * + (heatOfVap_density_term_box[vapBox] - + heatOfVap_density_term_box[liqBox]); + heatOfVap *= unit::K_TO_KJ_PER_MOL; } } +#endif } From 91bb6a97da086fca6a8eaeca37be5851ee71726b Mon Sep 17 00:00:00 2001 From: GregorySchwing Date: Fri, 10 Dec 2021 10:10:51 -0600 Subject: [PATCH 16/16] Added compressibility, enthalpy, and heatofVap to console --- src/ConsoleOutput.cpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/ConsoleOutput.cpp b/src/ConsoleOutput.cpp index df5a11240..960f2c04f 100644 --- a/src/ConsoleOutput.cpp +++ b/src/ConsoleOutput.cpp @@ -229,7 +229,11 @@ void ConsoleOutput::PrintStatistic(const uint box, const ulong step) const if(enablePressure) { printElement(var->pressure[box], elementWidth); - printElement((var->pressure[box]) * (var->volumeRef[box]) / (var->numByBox[box]) / (var->T_in_K) / (UNIT_CONST_H::unit::K_MOLECULE_PER_A3_TO_BAR), elementWidth); + printElement(var->compressibility[box], elementWidth); + printElement(var->enthalpy[box], elementWidth); + #if ENSEMBLE == GEMC + printElement(var->heatOfVap, elementWidth); + #endif } if(enableMol) { printElement(var->numByBox[box], elementWidth); @@ -340,6 +344,10 @@ void ConsoleOutput::PrintStatisticTitle() if(enablePressure) { printElement("PRESSURE", elementWidth); printElement("COMPRESSIBILITY", elementWidth); + printElement("ENTHALPY", elementWidth); + #if ENSEMBLE == GEMC + printElement("HEAT_VAP", elementWidth); + #endif } if(enableMol) {