diff --git a/plugins/openaccPTM/.gitignore b/plugins/openaccPTM/.gitignore new file mode 100644 index 0000000000..f11b086c1d --- /dev/null +++ b/plugins/openaccPTM/.gitignore @@ -0,0 +1,4 @@ +Makefile.conf +extra/ +.deps/ + diff --git a/plugins/openaccPTM/ACCParallelTaskManager.h b/plugins/openaccPTM/ACCParallelTaskManager.h new file mode 100644 index 0000000000..cc1480d7a3 --- /dev/null +++ b/plugins/openaccPTM/ACCParallelTaskManager.h @@ -0,0 +1,345 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2017-2023 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see . ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#ifndef __PLUMED_acc_ParallelTaskManager_h +#define __PLUMED_acc_ParallelTaskManager_h + +#include "plumed/core/ActionWithVector.h" +#include "plumed/core/ActionWithMatrix.h" +#include "plumed/core/ParallelTaskManager.h" +#include "plumed/tools/Communicator.h" +#include "plumed/tools/OpenMP.h" +#include "plumed/tools/View.h" +#include "plumed/tools/View2D.h" + +#include "plumed/tools/ColvarOutput.h" +#include "plumed/tools/OpenACC.h" + +namespace PLMD { + +template +class AccParallelTaskManager : public ParallelTaskManager { +public: + static void registerKeywords( Keywords& keys ) { + ParallelTaskManager::registerKeywords(keys); + } + AccParallelTaskManager(ActionWithVector* av): + PLMD::ParallelTaskManager(av) {} +/// This runs all the tasks + void runAllTasks(); +/// Apply the forces on the parallel object + void applyForces( std::vector& forcesForApply ); +}; + +struct ACCPTM { + template + using PTM=AccParallelTaskManager; +}; + +//use the __PLUMED_USE_OPENACC_TASKSMINE macro to debug the ptm ins a single file +//so that compiling witha a small modification will be faster (the ptm is included nearly everywhere) +#ifndef __PLUMED_USE_OPENACC_TASKSMINE +template +void runAllTasksACC(typename T::input_type actiondata, + ParallelActionsInput myinput, + std::vector& value_stash, + const std::vector & partialTaskList, + const unsigned nactive_tasks, + const std::size_t nderivatives_per_task, + const std::size_t workspace_size + ) { + auto myinput_acc = OpenACC::fromToDataHelper(myinput); + auto actiondata_acc = OpenACC::fromToDataHelper(actiondata); + + //template type is deduced + OpenACC::memoryManager vs{value_stash}; + auto value_stash_data = vs.devicePtr(); + + OpenACC::memoryManager ptl{partialTaskList}; + auto partialTaskList_data = ptl.devicePtr(); + + OpenACC::memoryManager buff{workspace_size*nactive_tasks}; + + auto buffer = buff.devicePtr(); + OpenACC::memoryManager dev(nderivatives_per_task*nactive_tasks); + auto derivatives = dev.devicePtr(); +#pragma acc parallel loop present(myinput, actiondata) \ + copyin(nactive_tasks, \ + nderivatives_per_task, \ + workspace_size)\ + deviceptr(derivatives, \ + partialTaskList_data, \ + value_stash_data, \ + buffer) \ + default(none) + for(unsigned i=0; i0)? + buffer+workspace_size*i + :nullptr ); + // Calculate the stuff in the loop for this action + T::performTask( task_index, actiondata, myinput, myout ); + } + vs.copyFromDevice(value_stash.data()); +} +#else +template +void runAllTasksACC(typename T::input_type actiondata, + ParallelActionsInput myinput, + std::vector& value_stash, + const std::vector & partialTaskList, + const unsigned nactive_tasks, + const std::size_t nderivatives_per_task, + const std::size_t workspace_size + ) ; +#endif //__PLUMED_USE_OPENACC_TASKSMINE + +template +void AccParallelTaskManager::runAllTasks() { + // Get the list of active tasks + std::vector & partialTaskList( this->action->getListOfActiveTasks( this->action ) ); + unsigned nactive_tasks=partialTaskList.size(); + // Get all the input data so we can broadcast it to the GPU + this->myinput.noderiv = true; + this->action->getInputData( this->input_buffer ); + this->myinput.dataSize = this->input_buffer.size(); + this->myinput.inputdata = this->input_buffer.data(); + // Transfer all the bookeeping information about the arguments + this->argumentsMap.setupArguments( this->action ); + this->myinput.setupArguments( this->argumentsMap ); + // Reset the values at the start of the task loop + std::size_t totalvals=this->getValueStashSize(); + if( this->value_stash.size()!=totalvals ) { + this->value_stash.resize(totalvals); + } + std::fill (this->value_stash.begin(),this->value_stash.end(), 0.0); + if (this->comm.Get_rank()== 0) {// no multigpu shenanigans until this works + runAllTasksACC( + this->actiondata, + this->myinput, + this->value_stash, + partialTaskList, + nactive_tasks, + this->nderivatives_per_task, + this->workspace_size + ); + } + this->comm.Bcast( this->value_stash.data(), this->value_stash.size(), 0); + + // MPI Gather everything + if( !this->action->runInSerial() ) { + this->comm.Sum( this->value_stash ); + } + +// And transfer the value to the output values + this->action->transferStashToValues( this->value_stash ); +} + +//use the __PLUMED_USE_OPENACC_FORCESMINE macro to debug the ptm ins a single file +//so that compiling witha a small modification will be faster (the ptm is included nearly everywhere) +#ifndef __PLUMED_USE_OPENACC_FORCESMINE +template +void applyForcesWithACC(PLMD::View forcesForApply, + typename T::input_type actiondata, + ParallelActionsInput myinput, + const std::vector& value_stash, + const std::vector & partialTaskList, + const unsigned nactive_tasks, + const std::size_t nderivatives_per_task, + const std::size_t workspace_size + ) { + auto myinput_acc = OpenACC::fromToDataHelper(myinput); + auto actiondata_acc = OpenACC::fromToDataHelper(actiondata); + + //template type is deduced + OpenACC::memoryManager vs{value_stash}; + auto value_stash_data = vs.devicePtr(); + + OpenACC::memoryManager ptl{partialTaskList}; + auto partialTaskList_data = ptl.devicePtr(); + + OpenACC::memoryManager ffa {forcesForApply}; + auto forcesForApply_data = ffa.devicePtr(); + const auto forcesForApply_size = ffa.size(); + const auto nind_per_scalar = ForceIndexHolder::indexesPerScalar(myinput); + //nscalars is >=ncomponents (see setupParallelTaskManager ) + const auto nind_per_task = nind_per_scalar*myinput.nscalars; + + OpenACC::memoryManager dev{nderivatives_per_task*nactive_tasks}; + auto derivatives = dev.devicePtr(); + OpenACC::memoryManager ind{nind_per_task*nactive_tasks}; + auto indices = ind.devicePtr(); + OpenACC::memoryManager vtmp{myinput.sizeOfFakeVals()*nactive_tasks}; + auto valstmp = vtmp.devicePtr(); + OpenACC::memoryManager buff{workspace_size*nactive_tasks}; + auto buffer = buff.devicePtr(); + +#define forces_indicesArg(taskID,scalarID) ForceIndexHolder::create(myinput, \ + indices + taskID*nind_per_task + scalarID*nind_per_scalar) +#define derivativeDrift(taskID,scalarID) taskID*nderivatives_per_task \ + + scalarID*myinput.ncomponents*myinput.nderivatives_per_scalar +#define stashDrift(taskID,scalarID) taskID*myinput.nscalars \ + + scalarID*myinput.ncomponents + +#pragma acc data present(myinput,actiondata) \ + copyin(nactive_tasks, \ + forcesForApply_size, \ + nderivatives_per_task, nind_per_task,nind_per_scalar, \ + workspace_size) \ + deviceptr(derivatives, \ + indices, \ + value_stash_data, \ + partialTaskList_data, \ + forcesForApply_data, \ + valstmp, \ + buffer) \ + default(none) + { +#pragma acc parallel loop + for(unsigned t=0; t0)?buffer+workspace_size*t:nullptr); + // Calculate the stuff in the loop for this action + T::performTask( task_index, actiondata, myinput, myout ); + // If this is a matrix this returns a number that isn't one as we have to loop over the columns + const std::size_t nvpt = T::getNumberOfValuesPerTask( task_index, actiondata ); +#pragma acc loop seq + for(unsigned vID=0; vID::gatherThreadSafeForces( myinput, + force_indices, + finput, + View(forcesForApply_data, + forcesForApply_size)); + } + } + +#pragma acc parallel loop + for(unsigned v=myinput.threadunsafe_forces_start; v +void applyForcesWithACC(PLMD::View forcesForApply, + typename T::input_type actiondata, + ParallelActionsInput myinput, + const std::vector& value_stash, + const std::vector & partialTaskList, + const unsigned nactive_tasks, + const std::size_t nderivatives_per_task, + const std::size_t workspace_size + ); +#endif //__PLUMED_USE_OPENACC_FORCESMINE + +template +void AccParallelTaskManager::applyForces( std::vector& forcesForApply ) { + // Get the list of active tasks + std::vector & partialTaskList= this->action->getListOfActiveTasks( this->action ) ; + unsigned nactive_tasks=partialTaskList.size(); + // Clear force buffer + forcesForApply.assign( forcesForApply.size(), 0.0 ); + //TODO: check if std::fill is faster (i get conflicting answers on the net) + //std::fill (forcesForApply.begin(),forcesForApply.end(), 0.0); + // Get all the input data so we can broadcast it to the GPU + this->myinput.noderiv = false; + // Retrieve the forces from the values + this->action->transferForcesToStash( this->value_stash ); + + std::fill (this->omp_forces[0].begin(),this->omp_forces[0].end(), 0.0); + if (this->comm.Get_rank() == 0) { + applyForcesWithACC( + PLMD::View { forcesForApply.data(), forcesForApply.size() }, + this->actiondata, + this->myinput, + this->value_stash, + partialTaskList, + nactive_tasks, + this->nderivatives_per_task, + this->workspace_size + ); + } + // MPI Gather everything (this must be extended to the gpu thing, after makning it mpi-aware) + if( !this->action->runInSerial() ) { + this->comm.Sum( forcesForApply ); + } +} + +} // namespace PLMD +#endif // __PLUMED_acc_ParallelTaskManager_h + + diff --git a/plugins/openaccPTM/Angle.cpp b/plugins/openaccPTM/Angle.cpp new file mode 100644 index 0000000000..4779e6b121 --- /dev/null +++ b/plugins/openaccPTM/Angle.cpp @@ -0,0 +1,30 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2011-2023 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see . ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#include "plumed/core/ActionRegister.h" +#include "plumed/colvar/MultiColvarTemplate.h" +#include "plumed/colvar/Angle.h" + +#include "ACCParallelTaskManager.h" + +typedef PLMD::colvar::MultiColvarTemplate AngleMultiAcc; + +PLUMED_REGISTER_ACTION(AngleMultiAcc,"ANGLE_VECTORACC") diff --git a/plugins/openaccPTM/DihedralCorrelation.cpp b/plugins/openaccPTM/DihedralCorrelation.cpp new file mode 100644 index 0000000000..f77585c304 --- /dev/null +++ b/plugins/openaccPTM/DihedralCorrelation.cpp @@ -0,0 +1,30 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2011-2023 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see . ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#include "plumed/core/ActionRegister.h" +#include "plumed/colvar/MultiColvarTemplate.h" +#include "plumed/colvar/DihedralCorrelation.h" + +#include "ACCParallelTaskManager.h" + +typedef PLMD::colvar::MultiColvarTemplate DihedralCorrelationMultiAcc; + +PLUMED_REGISTER_ACTION(DihedralCorrelationMultiAcc,"DIHEDRALCORRELATION_VECTORACC") diff --git a/plugins/openaccPTM/Dipole.cpp b/plugins/openaccPTM/Dipole.cpp new file mode 100644 index 0000000000..af5801e1a0 --- /dev/null +++ b/plugins/openaccPTM/Dipole.cpp @@ -0,0 +1,30 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2011-2023 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see . ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#include "plumed/core/ActionRegister.h" +#include "plumed/colvar/MultiColvarTemplate.h" +#include "plumed/colvar/Dipole.h" + +#include "ACCParallelTaskManager.h" + +typedef PLMD::colvar::MultiColvarTemplate DipoleMultiAcc; + +PLUMED_REGISTER_ACTION(DipoleMultiAcc,"DIPOLE_VECTORACC") diff --git a/plugins/openaccPTM/Distance.cpp b/plugins/openaccPTM/Distance.cpp new file mode 100644 index 0000000000..03d4c21dfc --- /dev/null +++ b/plugins/openaccPTM/Distance.cpp @@ -0,0 +1,30 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2011-2023 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see . ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#include "plumed/core/ActionRegister.h" +#include "plumed/colvar/MultiColvarTemplate.h" +#include "plumed/colvar/Distance.h" + +#include "ACCParallelTaskManager.h" + +typedef PLMD::colvar::MultiColvarTemplate DistanceMultiAcc; + +PLUMED_REGISTER_ACTION(DistanceMultiAcc,"DISTANCE_VECTORACC") diff --git a/plugins/openaccPTM/Makefile b/plugins/openaccPTM/Makefile new file mode 100644 index 0000000000..8d2e9391da --- /dev/null +++ b/plugins/openaccPTM/Makefile @@ -0,0 +1,67 @@ +include ./Makefile.conf + +#Dependency tracking based on https://make.mad-scientist.net/papers/advanced-auto-dependency-generation/#tldr +#this assumes gcc +DEPDIR := .deps +DEPFLAGS = -MT $@ -MMD -MP -MF $(DEPDIR)/$*.d +#tested with nvcc with :"Build cuda_11.7.r11.7/compiler.31442593_0" +#-dc adds relocatable device code +#-dlto Perform link-time optimization of device code. +OBJS = Angle.o DihedralCorrelation.o Dipole.o Distance.o Plane.o Position.o Torsion.o extra/colvar/ColvarInput.o extra/tools/Pbc.o extra/tools/Torsion.o extra/tools/Angle.o + +ADDCPPFLAGS=$(PLUMED_INCLUDE) -D__PLUMED_HAS_OPENACC=1 -acc -acc="gpu" +ADDCLDFLAGS=-acc -static-nvidia + +ifeq ($(SOEXT),dylib) + SONAME_OPTION:=-Xlinker -install_name +else + SONAME_OPTION:=-Xlinker -soname +endif + +PLUGIN=plumedOpenACC.$(SOEXT) + +NVCXX ?= nvc++ +NVLDSHARED = $(LDSHARED:$(CXX)=$(NVCXX)) +all: $(PLUGIN) + +%.o: %.cpp $(DEPDIR)/%.d | $(DEPDIR) + @echo Compiling object $@ + $(NVCXX) -c $(DEPFLAGS) $(CPPFLAGS) $(ADDCPPFLAGS) $(CXXFLAGS) $< -o $@ + +extra/colvar/ColvarInput.o: plumed/colvar/ColvarInput.cpp $(DEPDIR)/extra/colvar/ColvarInput.d | $(DEPDIR)/extra/colvar extra/colvar + @echo Compiling object $@ + $(NVCXX) -c $(DEPFLAGS) $(CPPFLAGS) $(ADDCPPFLAGS) $(CXXFLAGS) $< -o $@ + +extra/tools/Pbc.o: plumed/tools/Pbc.cpp $(DEPDIR)/extra/tools/Pbc.d | $(DEPDIR)/extra/tools extra/tools + @echo Compiling object $@ + $(NVCXX) -c $(DEPFLAGS) $(CPPFLAGS) $(ADDCPPFLAGS) $(CXXFLAGS) $< -o $@ + +extra/tools/Torsion.o: plumed/tools/Torsion.cpp $(DEPDIR)/extra/tools/Torsion.d | $(DEPDIR)/extra/tools extra/tools + @echo Compiling object $@ + $(NVCXX) -c $(DEPFLAGS) $(CPPFLAGS) $(ADDCPPFLAGS) $(CXXFLAGS) $< -o $@ + +extra/tools/Angle.o: plumed/tools/Angle.cpp $(DEPDIR)/extra/tools/Angle.d | $(DEPDIR)/extra/tools extra/tools + @echo Compiling object $@ + $(NVCXX) -c $(DEPFLAGS) $(CPPFLAGS) $(ADDCPPFLAGS) $(CXXFLAGS) $< -o $@ + +$(DEPDIR): ; @mkdir -p $@ +$(DEPDIR)/extra/colvar: ; @mkdir -p $@ +$(DEPDIR)/extra/tools: ; @mkdir -p $@ +extra/colvar: ; mkdir -p $@ +extra/tools: ; mkdir -p $@ + +DEPFILES := $(OBJS:%.o=$(DEPDIR)/%.d) +$(DEPFILES): +include $(wildcard $(DEPFILES)) + +$(PLUGIN): $(OBJS) + @echo Linking $@ + $(NVLDSHARED) $(ADDCLDFLAGS) $(SONAME_OPTION),"$(notdir $@)" $(DYNAMIC_LIBS) $(PLUMED_KERNEL) -o $@ $^ + +clean: + @rm -fv $(OBJS) $(PLUGIN) + @rm -rf .deps + +check: all + $(MAKE) -C regtest + $(MAKE) -C regtest checkfail diff --git a/plugins/openaccPTM/Plane.cpp b/plugins/openaccPTM/Plane.cpp new file mode 100644 index 0000000000..a771f4980a --- /dev/null +++ b/plugins/openaccPTM/Plane.cpp @@ -0,0 +1,30 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2011-2023 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see . ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#include "plumed/core/ActionRegister.h" +#include "plumed/colvar/MultiColvarTemplate.h" +#include "plumed/colvar/Plane.h" + +#include "ACCParallelTaskManager.h" + +typedef PLMD::colvar::MultiColvarTemplate PlaneMultiAcc; + +PLUMED_REGISTER_ACTION(PlaneMultiAcc,"PLANE_VECTORACC") diff --git a/plugins/openaccPTM/Position.cpp b/plugins/openaccPTM/Position.cpp new file mode 100644 index 0000000000..41bd8e6e44 --- /dev/null +++ b/plugins/openaccPTM/Position.cpp @@ -0,0 +1,30 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2011-2023 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see . ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#include "plumed/core/ActionRegister.h" +#include "plumed/colvar/MultiColvarTemplate.h" +#include "plumed/colvar/Position.h" + +#include "ACCParallelTaskManager.h" + +typedef PLMD::colvar::MultiColvarTemplate PositionMultiAcc; + +PLUMED_REGISTER_ACTION(PositionMultiAcc,"POSITION_VECTORACC") diff --git a/plugins/openaccPTM/Torsion.cpp b/plugins/openaccPTM/Torsion.cpp new file mode 100644 index 0000000000..709680057c --- /dev/null +++ b/plugins/openaccPTM/Torsion.cpp @@ -0,0 +1,30 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2011-2023 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see . ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#include "plumed/core/ActionRegister.h" +#include "plumed/colvar/MultiColvarTemplate.h" +#include "plumed/colvar/Torsion.h" + +#include "ACCParallelTaskManager.h" + +typedef PLMD::colvar::MultiColvarTemplate TorsionMultiAcc; + +PLUMED_REGISTER_ACTION(TorsionMultiAcc,"TORSION_VECTORACC") diff --git a/plugins/openaccPTM/configure.sh b/plugins/openaccPTM/configure.sh new file mode 100755 index 0000000000..f6651959b1 --- /dev/null +++ b/plugins/openaccPTM/configure.sh @@ -0,0 +1,30 @@ +#! /usr/bin/env bash + +if [[ -z $PLUMED_KERNEL ]]; then + echo "$(basename $0) can work only if \"PLUMED_KERNEL\" is defined" + echo "either via module load or sourceme.sh" +fi + +{ + plumed config makefile_conf + echo "PLUMED_INCLUDE=-I$(plumed info --include-dir)" + echo "PLUMED_KERNEL=-L${PLUMED_KERNEL}" +} >Make.tmp + +#pendantic adds a unuseful FOR EACH line with +#"" warning: style of line directive is a GCC extension" +{ + # grep CXXFLAGS Make.tmp \ + # | sed -e 's/-f/-Xcompiler -f/g' \ + # -e 's/-pedantic//g' \ + # -e 's/-W/-Xcompiler -W/g' + # grep -eDYNAMIC_LIBS -eLDFLAGS Make.tmp \ + # | sed -e 's/-rdynamic/-Xcompiler -rdynamic/g' \ + # -e 's/-Wl,/-Xlinker /g' \ + # -e 's/-f/-Xcompiler -f/g' + # #prints the rest of the file + # grep -eDYNAMIC_LIBS -eLDFLAGS -eCXXFLAGS Make.tmp -v + sed -e 's/-fno-gnu-unique//g' Make.tmp +} >Makefile.conf + +rm Make.tmp diff --git a/plugins/openaccPTM/plumed b/plugins/openaccPTM/plumed new file mode 120000 index 0000000000..929cb3dc9b --- /dev/null +++ b/plugins/openaccPTM/plumed @@ -0,0 +1 @@ +../../src \ No newline at end of file diff --git a/src/colvar/Angle.cpp b/src/colvar/Angle.cpp index b0e0866725..473cd5addd 100644 --- a/src/colvar/Angle.cpp +++ b/src/colvar/Angle.cpp @@ -19,10 +19,7 @@ You should have received a copy of the GNU Lesser General Public License along with plumed. If not, see . +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#ifdef __PLUMED_HAS_OPENACC -#define __PLUMED_USE_OPENACC 1 -#endif //__PLUMED_HAS_OPENACC -#include "Colvar.h" +#include "Angle.h" #include "ColvarShortcut.h" #include "MultiColvarTemplate.h" #include "core/ActionRegister.h" @@ -100,108 +97,12 @@ PRINT ARG=a4,a5,a6 FILE=colvar */ //+ENDPLUMEDOC -class Angle : public Colvar { - bool pbc; - std::vector value; - std::vector derivs; -public: - explicit Angle(const ActionOptions&); -// active methods: - void calculate() override; - static void registerKeywords( Keywords& keys ); - static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); - static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ); -}; - typedef ColvarShortcut AngleShortcut; PLUMED_REGISTER_ACTION(AngleShortcut,"ANGLE") PLUMED_REGISTER_ACTION(Angle,"ANGLE_SCALAR") typedef MultiColvarTemplate AngleMulti; PLUMED_REGISTER_ACTION(AngleMulti,"ANGLE_VECTOR") -void Angle::registerKeywords( Keywords& keys ) { - Colvar::registerKeywords(keys); - keys.setDisplayName("ANGLE"); - keys.add("atoms","ATOMS","the list of atoms involved in this collective variable (either 3 or 4 atoms)"); - keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log"); - keys.setValueDescription("scalar/vector","the ANGLE involving these atoms"); - keys.reset_style("NUMERICAL_DERIVATIVES","hidden"); -} - -void Angle::parseAtomList( const int& num, std::vector& atoms, ActionAtomistic* aa ) { - aa->parseAtomList("ATOMS",num,atoms); - if(atoms.size()==3) { - aa->log.printf(" between atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial()); - atoms.resize(4); - atoms[3]=atoms[2]; - atoms[2]=atoms[1]; - } else if(atoms.size()==4) { - aa->log.printf(" between lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial()); - } else if( num<0 || atoms.size()>0 ) { - aa->error("Number of specified atoms should be either 3 or 4"); - } -} - -unsigned Angle::getModeAndSetupValues( ActionWithValue* av ) { - av->addValueWithDerivatives(); - av->setNotPeriodic(); - return 0; -} - -Angle::Angle(const ActionOptions&ao): - PLUMED_COLVAR_INIT(ao), - pbc(true), - value(1) { - std::vector atoms; - parseAtomList( -1, atoms, this ); - bool nopbc=!pbc; - parseFlag("NOPBC",nopbc); - pbc=!nopbc; - - if(pbc) { - log.printf(" using periodic boundary conditions\n"); - } else { - log.printf(" without periodic boundary conditions\n"); - } - - addValueWithDerivatives(); - setNotPeriodic(); - requestAtoms(atoms); - checkRead(); -} - -// calculator -void Angle::calculate() { - - if(pbc) { - makeWhole(); - } - ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this); - calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout ); - setValue( value[0] ); - for(unsigned i=0; i. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#ifndef __PLUMED_colvar_Angle_h +#define __PLUMED_colvar_Angle_h +#include "Colvar.h" +#include "ColvarInput.h" +#include "tools/Angle.h" + +namespace PLMD { +namespace colvar { + +class Angle : public Colvar { + bool pbc; + std::vector value; + std::vector derivs; +public: + explicit Angle(const ActionOptions&); +// active methods: + void calculate() override; + static void registerKeywords( Keywords& keys ); + static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); + static unsigned getModeAndSetupValues( ActionWithValue* av ); + static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ); +}; + +void Angle::registerKeywords( Keywords& keys ) { + Colvar::registerKeywords(keys); + keys.setDisplayName("ANGLE"); + keys.add("atoms","ATOMS","the list of atoms involved in this collective variable (either 3 or 4 atoms)"); + keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log"); + keys.setValueDescription("scalar/vector","the ANGLE involving these atoms"); + keys.reset_style("NUMERICAL_DERIVATIVES","hidden"); +} + +void Angle::parseAtomList( const int& num, std::vector& atoms, ActionAtomistic* aa ) { + aa->parseAtomList("ATOMS",num,atoms); + if(atoms.size()==3) { + aa->log.printf(" between atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial()); + atoms.resize(4); + atoms[3]=atoms[2]; + atoms[2]=atoms[1]; + } else if(atoms.size()==4) { + aa->log.printf(" between lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial()); + } else if( num<0 || atoms.size()>0 ) { + aa->error("Number of specified atoms should be either 3 or 4"); + } +} + +unsigned Angle::getModeAndSetupValues( ActionWithValue* av ) { + av->addValueWithDerivatives(); + av->setNotPeriodic(); + return 0; +} + +Angle::Angle(const ActionOptions&ao): + PLUMED_COLVAR_INIT(ao), + pbc(true), + value(1) { + std::vector atoms; + parseAtomList( -1, atoms, this ); + bool nopbc=!pbc; + parseFlag("NOPBC",nopbc); + pbc=!nopbc; + + if(pbc) { + log.printf(" using periodic boundary conditions\n"); + } else { + log.printf(" without periodic boundary conditions\n"); + } + + addValueWithDerivatives(); + setNotPeriodic(); + requestAtoms(atoms); + checkRead(); +} + +// calculator +void Angle::calculate() { + + if(pbc) { + makeWhole(); + } + ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this); + calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout ); + setValue( value[0] ); + for(unsigned i=0; i::registerKeywords(Keywords& keys ) { } keys.addActionNameSuffix("_SCALAR"); keys.addActionNameSuffix("_VECTOR"); + //GPU related settings + keys.addFlag("USEGPU",false,"run this calculation on the GPU"); + keys.addLinkInDocForFlag("USEGPU","gpu.md"); + keys.addActionNameSuffix("_VECTORACC"); } template @@ -64,8 +68,12 @@ ColvarShortcut::ColvarShortcut(const ActionOptions&ao): if( keywords.style( key, "atoms" ) ) { std::string inpt; parseNumbered( key, 1, inpt ); + bool usegpuFLAG=false; + parseFlag("USEGPU",usegpuFLAG); if( inpt.length()>0 ) { - readInputLine( getShortcutLabel() + ": " + getName() + "_VECTOR " + key + "1=" + inpt + " " + convertInputLineToString() ); + readInputLine( getShortcutLabel() + ": " + + getName() + "_VECTOR" + (usegpuFLAG ? "ACC ":" ") + + key + "1=" + inpt + " " + convertInputLineToString() ); scalar=false; break; } diff --git a/src/colvar/DihedralCorrelation.cpp b/src/colvar/DihedralCorrelation.cpp index b5a93ff599..65cb5598ab 100644 --- a/src/colvar/DihedralCorrelation.cpp +++ b/src/colvar/DihedralCorrelation.cpp @@ -19,18 +19,12 @@ You should have received a copy of the GNU Lesser General Public License along with plumed. If not, see . +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#ifdef __PLUMED_HAS_OPENACC -#define __PLUMED_USE_OPENACC 1 -#endif //__PLUMED_HAS_OPENACC -#include "Colvar.h" +#include "DihedralCorrelation.h" #include "ColvarShortcut.h" #include "MultiColvarTemplate.h" #include "tools/Torsion.h" #include "core/ActionRegister.h" -#include -#include - namespace PLMD { namespace colvar { @@ -80,19 +74,6 @@ PRINT ARG=d FILE=colvar */ //+ENDPLUMEDOC -class DihedralCorrelation : public Colvar { -private: - bool pbc; - std::vector value; - std::vector derivs; -public: - static void registerKeywords( Keywords& keys ); - explicit DihedralCorrelation(const ActionOptions&); - static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); - void calculate() override; - static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ); -}; typedef ColvarShortcut DihedralCorrelationShortcut; PLUMED_REGISTER_ACTION(DihedralCorrelationShortcut,"DIHEDRAL_CORRELATION") @@ -100,112 +81,5 @@ PLUMED_REGISTER_ACTION(DihedralCorrelation,"DIHEDRAL_CORRELATION_SCALAR") typedef MultiColvarTemplate DihedralCorrelationMulti; PLUMED_REGISTER_ACTION(DihedralCorrelationMulti,"DIHEDRAL_CORRELATION_VECTOR") -void DihedralCorrelation::registerKeywords( Keywords& keys ) { - Colvar::registerKeywords( keys ); - keys.setDisplayName("DIHEDRAL_CORRELATION"); - keys.add("atoms","ATOMS","the set of 8 atoms that are being used to calculate this quantity"); - keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log"); - keys.setValueDescription("scalar/vector","the DIHEDRAL_CORRELATION for these atoms"); - keys.reset_style("NUMERICAL_DERIVATIVES","hidden"); -} - -DihedralCorrelation::DihedralCorrelation(const ActionOptions&ao): - PLUMED_COLVAR_INIT(ao), - pbc(true), - value(1) { - std::vector atoms; - parseAtomList(-1,atoms,this); - if( atoms.size()!=8 ) { - error("Number of specified atoms should be 8"); - } - - bool nopbc=!pbc; - parseFlag("NOPBC",nopbc); - pbc=!nopbc; - - if(pbc) { - log.printf(" using periodic boundary conditions\n"); - } else { - log.printf(" without periodic boundary conditions\n"); - } - addValueWithDerivatives(); - setNotPeriodic(); -} - -void DihedralCorrelation::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { - aa->parseAtomList("ATOMS",num,t); - if( num<0 && t.size()!=8 ) { - aa->error("Number of specified atoms should be 8"); - } - if( t.size()==8 ) { - aa->log.printf(" correlation between dihedral angle for atoms %d %d %d %d and atoms %d %d %d %d\n", - t[0].serial(),t[1].serial(),t[2].serial(),t[3].serial(),t[4].serial(),t[5].serial(),t[6].serial(),t[7].serial()); - } -} - -unsigned DihedralCorrelation::getModeAndSetupValues( ActionWithValue* av ) { - av->addValueWithDerivatives(); - av->setNotPeriodic(); - return 0; -} - -void DihedralCorrelation::calculate() { - - if(pbc) { - makeWhole(); - } - ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this); - calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout ); - setValue( value[0] ); - for(unsigned i=0; i. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#ifndef __PLUMED_colvar_DihedralCorrelation_h +#define __PLUMED_colvar_DihedralCorrelation_h +#include "Colvar.h" +#include "ColvarInput.h" +#include "tools/Torsion.h" + +#include +#include + +namespace PLMD { +namespace colvar { + +class DihedralCorrelation : public Colvar { +private: + bool pbc; + std::vector value; + std::vector derivs; +public: + static void registerKeywords( Keywords& keys ); + explicit DihedralCorrelation(const ActionOptions&); + static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); + static unsigned getModeAndSetupValues( ActionWithValue* av ); + void calculate() override; + static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ); +}; + +void DihedralCorrelation::registerKeywords( Keywords& keys ) { + Colvar::registerKeywords( keys ); + keys.setDisplayName("DIHEDRAL_CORRELATION"); + keys.add("atoms","ATOMS","the set of 8 atoms that are being used to calculate this quantity"); + keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log"); + keys.setValueDescription("scalar/vector","the DIHEDRAL_CORRELATION for these atoms"); + keys.reset_style("NUMERICAL_DERIVATIVES","hidden"); +} + +DihedralCorrelation::DihedralCorrelation(const ActionOptions&ao): + PLUMED_COLVAR_INIT(ao), + pbc(true), + value(1) { + std::vector atoms; + parseAtomList(-1,atoms,this); + if( atoms.size()!=8 ) { + error("Number of specified atoms should be 8"); + } + + bool nopbc=!pbc; + parseFlag("NOPBC",nopbc); + pbc=!nopbc; + + if(pbc) { + log.printf(" using periodic boundary conditions\n"); + } else { + log.printf(" without periodic boundary conditions\n"); + } + addValueWithDerivatives(); + setNotPeriodic(); +} + +void DihedralCorrelation::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { + aa->parseAtomList("ATOMS",num,t); + if( num<0 && t.size()!=8 ) { + aa->error("Number of specified atoms should be 8"); + } + if( t.size()==8 ) { + aa->log.printf(" correlation between dihedral angle for atoms %d %d %d %d and atoms %d %d %d %d\n", + t[0].serial(),t[1].serial(),t[2].serial(),t[3].serial(),t[4].serial(),t[5].serial(),t[6].serial(),t[7].serial()); + } +} + +unsigned DihedralCorrelation::getModeAndSetupValues( ActionWithValue* av ) { + av->addValueWithDerivatives(); + av->setNotPeriodic(); + return 0; +} + +void DihedralCorrelation::calculate() { + + if(pbc) { + makeWhole(); + } + ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this); + calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout ); + setValue( value[0] ); + for(unsigned i=0; i. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#ifdef __PLUMED_HAS_OPENACC -#define __PLUMED_USE_OPENACC 1 -#endif //__PLUMED_HAS_OPENACC -#include "Colvar.h" +#include "Dipole.h" #include "ColvarShortcut.h" #include "MultiColvarTemplate.h" #include "core/ActionRegister.h" @@ -91,165 +88,11 @@ PRINT FILE=output STRIDE=5 ARG=d */ //+ENDPLUMEDOC -class Dipole : public Colvar { - std::vector ga_lista; - bool components; - bool nopbc; - std::vector value; - std::vector derivs; - Value* valuex=nullptr; - Value* valuey=nullptr; - Value* valuez=nullptr; -public: - explicit Dipole(const ActionOptions&); - static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); - void calculate() override; - static void registerKeywords(Keywords& keys); - static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ); -}; - typedef ColvarShortcut DipoleShortcut; PLUMED_REGISTER_ACTION(DipoleShortcut,"DIPOLE") PLUMED_REGISTER_ACTION(Dipole,"DIPOLE_SCALAR") typedef MultiColvarTemplate DipoleMulti; PLUMED_REGISTER_ACTION(DipoleMulti,"DIPOLE_VECTOR") -void Dipole::registerKeywords(Keywords& keys) { - Colvar::registerKeywords(keys); - keys.setDisplayName("DIPOLE"); - keys.add("atoms","GROUP","the group of atoms we are calculating the dipole moment for"); - keys.addFlag("COMPONENTS",false,"calculate the x, y and z components of the dipole separately and store them as label.x, label.y and label.z"); - keys.addOutputComponent("x","COMPONENTS","scalar/vector","the x-component of the dipole"); - keys.addOutputComponent("y","COMPONENTS","scalar/vector","the y-component of the dipole"); - keys.addOutputComponent("z","COMPONENTS","scalar/vector","the z-component of the dipole"); - keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log"); - keys.setValueDescription("scalar/vector","the DIPOLE for these atoms"); - keys.reset_style("NUMERICAL_DERIVATIVES","hidden"); -} - -Dipole::Dipole(const ActionOptions&ao): - PLUMED_COLVAR_INIT(ao), - components(false), - value(1) { - parseAtomList(-1,ga_lista,this); - components=(getModeAndSetupValues(this)==1); - if( components ) { - value.resize(3); - valuex=getPntrToComponent("x"); - valuey=getPntrToComponent("y"); - valuez=getPntrToComponent("z"); - } else { - derivs.resize(1,ga_lista.size()); - } - parseFlag("NOPBC",nopbc); - checkRead(); - - if(nopbc) { - log.printf(" without periodic boundary conditions\n"); - } else { - log.printf(" using periodic boundary conditions\n"); - } - - requestAtoms(ga_lista); -} - -void Dipole::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { - aa->parseAtomList("GROUP",num,t); - if( t.size()>0 ) { - aa->log.printf(" of %u atoms\n",static_cast(t.size())); - for(unsigned int i=0; ilog.printf(" %d", t[i].serial()); - } - aa->log.printf(" \n"); - } -} - -unsigned Dipole::getModeAndSetupValues( ActionWithValue* av ) { - bool c; - av->parseFlag("COMPONENTS",c); - if( c ) { - av->addComponentWithDerivatives("x"); - av->componentIsNotPeriodic("x"); - av->addComponentWithDerivatives("y"); - av->componentIsNotPeriodic("y"); - av->addComponentWithDerivatives("z"); - av->componentIsNotPeriodic("z"); - return 1; - } - av->addValueWithDerivatives(); - av->setNotPeriodic(); - return 0; -} - -// calculator -void Dipole::calculate() { - if( !chargesWereSet ) { - error("charges were not set by MD code"); - } - - if(!nopbc) { - makeWhole(); - } - unsigned N=getNumberOfAtoms(); - - if(!components) { - ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) ); - calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout ); - for(unsigned i=0; iset(value[0]); - valuey->set(value[1]); - valuez->set(value[2]); - } -} - -void Dipole::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) { - unsigned N=cvin.pos.size(); - double ctot=0.; - for(unsigned i=0; i. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#ifndef __PLUMED_colvar_Dipole_h +#define __PLUMED_colvar_Dipole_h +#include "Colvar.h" +#include "ColvarInput.h" + +namespace PLMD { +namespace colvar { + +class Dipole : public Colvar { + std::vector ga_lista; + bool components; + bool nopbc; + std::vector value; + std::vector derivs; + Value* valuex=nullptr; + Value* valuey=nullptr; + Value* valuez=nullptr; +public: + explicit Dipole(const ActionOptions&); + static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); + static unsigned getModeAndSetupValues( ActionWithValue* av ); + void calculate() override; + static void registerKeywords(Keywords& keys); + static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ); +}; + +void Dipole::registerKeywords(Keywords& keys) { + Colvar::registerKeywords(keys); + keys.setDisplayName("DIPOLE"); + keys.add("atoms","GROUP","the group of atoms we are calculating the dipole moment for"); + keys.addFlag("COMPONENTS",false,"calculate the x, y and z components of the dipole separately and store them as label.x, label.y and label.z"); + keys.addOutputComponent("x","COMPONENTS","scalar/vector","the x-component of the dipole"); + keys.addOutputComponent("y","COMPONENTS","scalar/vector","the y-component of the dipole"); + keys.addOutputComponent("z","COMPONENTS","scalar/vector","the z-component of the dipole"); + keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log"); + keys.setValueDescription("scalar/vector","the DIPOLE for these atoms"); + keys.reset_style("NUMERICAL_DERIVATIVES","hidden"); +} + +Dipole::Dipole(const ActionOptions&ao): + PLUMED_COLVAR_INIT(ao), + components(false), + value(1) { + parseAtomList(-1,ga_lista,this); + components=(getModeAndSetupValues(this)==1); + if( components ) { + value.resize(3); + valuex=getPntrToComponent("x"); + valuey=getPntrToComponent("y"); + valuez=getPntrToComponent("z"); + } else { + derivs.resize(1,ga_lista.size()); + } + parseFlag("NOPBC",nopbc); + checkRead(); + + if(nopbc) { + log.printf(" without periodic boundary conditions\n"); + } else { + log.printf(" using periodic boundary conditions\n"); + } + + requestAtoms(ga_lista); +} + +void Dipole::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { + aa->parseAtomList("GROUP",num,t); + if( t.size()>0 ) { + aa->log.printf(" of %u atoms\n",static_cast(t.size())); + for(unsigned int i=0; ilog.printf(" %d", t[i].serial()); + } + aa->log.printf(" \n"); + } +} + +unsigned Dipole::getModeAndSetupValues( ActionWithValue* av ) { + bool c; + av->parseFlag("COMPONENTS",c); + if( c ) { + av->addComponentWithDerivatives("x"); + av->componentIsNotPeriodic("x"); + av->addComponentWithDerivatives("y"); + av->componentIsNotPeriodic("y"); + av->addComponentWithDerivatives("z"); + av->componentIsNotPeriodic("z"); + return 1; + } + av->addValueWithDerivatives(); + av->setNotPeriodic(); + return 0; +} + +// calculator +void Dipole::calculate() { + if( !chargesWereSet ) { + error("charges were not set by MD code"); + } + + if(!nopbc) { + makeWhole(); + } + unsigned N=getNumberOfAtoms(); + + if(!components) { + ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) ); + calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout ); + for(unsigned i=0; iset(value[0]); + valuey->set(value[1]); + valuez->set(value[2]); + } +} + +void Dipole::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) { + unsigned N=cvin.pos.size(); + double ctot=0.; + for(unsigned i=0; i. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -//this is temporary: -#ifdef __PLUMED_HAS_OPENACC -#define __PLUMED_USE_OPENACC 1 -#endif //__PLUMED_HAS_OPENACC -#include "Colvar.h" + +#include "Distance.h" #include "ColvarShortcut.h" -#include "MultiColvarTemplate.h" #include "core/ActionRegister.h" -#include "tools/Pbc.h" namespace PLMD { namespace colvar { @@ -142,221 +137,11 @@ with domain (-0.5,+0.5). */ //+ENDPLUMEDOC -class Distance : public Colvar { - bool components; - bool scaled_components; - bool pbc; - - std::vector value; - std::vector derivs; -public: - static void registerKeywords( Keywords& keys ); - explicit Distance(const ActionOptions&); - static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); -// active methods: - void calculate() override; - static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ); -}; - typedef ColvarShortcut DistanceShortcut; PLUMED_REGISTER_ACTION(DistanceShortcut,"DISTANCE") PLUMED_REGISTER_ACTION(Distance,"DISTANCE_SCALAR") typedef MultiColvarTemplate DistanceMulti; PLUMED_REGISTER_ACTION(DistanceMulti,"DISTANCE_VECTOR") -void Distance::registerKeywords( Keywords& keys ) { - Colvar::registerKeywords( keys ); - keys.setDisplayName("DISTANCE"); - constexpr auto scalarOrVector = Keywords::componentType::scalar | Keywords::componentType::vector; - keys.add("atoms","ATOMS","the pair of atom that we are calculating the distance between"); - keys.addFlag("COMPONENTS",false,"calculate the x, y and z components of the distance separately and store them as label.x, label.y and label.z"); - keys.addFlag("SCALED_COMPONENTS",false,"calculate the a, b and c scaled components of the distance separately and store them as label.a, label.b and label.c"); - keys.addOutputComponent("x","COMPONENTS",scalarOrVector,"the x-component of the vector connecting the two atoms"); - keys.addOutputComponent("y","COMPONENTS",scalarOrVector,"the y-component of the vector connecting the two atoms"); - keys.addOutputComponent("z","COMPONENTS",scalarOrVector,"the z-component of the vector connecting the two atoms"); - keys.addOutputComponent("a","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the first lattice vector of the vector connecting the two atoms"); - keys.addOutputComponent("b","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the second lattice vector of the vector connecting the two atoms"); - keys.addOutputComponent("c","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the third lattice vector of the vector connecting the two atoms"); - keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log"); - keys.setValueDescription(scalarOrVector,"the DISTANCE between this pair of atoms"); - keys.addDOI("10.1007/978-1-4939-9608-7_21"); - keys.reset_style("NUMERICAL_DERIVATIVES","hidden"); -} - -Distance::Distance(const ActionOptions&ao): - PLUMED_COLVAR_INIT(ao), - components(false), - scaled_components(false), - pbc(true), - value(1) { - std::vector atoms; - parseAtomList(-1,atoms,this); - if(atoms.size()!=2) { - error("Number of specified atoms should be 2"); - } - - bool nopbc=!pbc; - parseFlag("NOPBC",nopbc); - pbc=!nopbc; - - if(pbc) { - log.printf(" using periodic boundary conditions\n"); - } else { - log.printf(" without periodic boundary conditions\n"); - } - - unsigned mode = getModeAndSetupValues( this ); - if(mode==1) { - components=true; - } else if(mode==2) { - scaled_components=true; - } - if( components || scaled_components ) { - value.resize(3); - } - requestAtoms(atoms); -} - -void Distance::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { - aa->parseAtomList("ATOMS",num,t); - if( t.size()==2 ) { - aa->log.printf(" between atoms %d %d\n",t[0].serial(),t[1].serial()); - } -} - -unsigned Distance::getModeAndSetupValues( ActionWithValue* av ) { - bool c; - av->parseFlag("COMPONENTS",c); - bool sc; - av->parseFlag("SCALED_COMPONENTS",sc); - if( c && sc ) { - av->error("COMPONENTS and SCALED_COMPONENTS are not compatible"); - } - - if(c) { - av->addComponentWithDerivatives("x"); - av->componentIsNotPeriodic("x"); - av->addComponentWithDerivatives("y"); - av->componentIsNotPeriodic("y"); - av->addComponentWithDerivatives("z"); - av->componentIsNotPeriodic("z"); - av->log<<" WARNING: components will not have the proper periodicity - see manual\n"; - return 1; - } else if(sc) { - av->addComponentWithDerivatives("a"); - av->componentIsPeriodic("a","-0.5","+0.5"); - av->addComponentWithDerivatives("b"); - av->componentIsPeriodic("b","-0.5","+0.5"); - av->addComponentWithDerivatives("c"); - av->componentIsPeriodic("c","-0.5","+0.5"); - return 2; - } - av->addValueWithDerivatives(); - av->setNotPeriodic(); - return 0; -} - -// calculator -void Distance::calculate() { - - if(pbc) { - makeWhole(); - } - - if( components ) { - ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) ); - calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), cvout ); - Value* valuex=getPntrToComponent("x"); - Value* valuey=getPntrToComponent("y"); - Value* valuez=getPntrToComponent("z"); - - for(unsigned i=0; i<2; ++i) { - setAtomsDerivatives(valuex,i,cvout.getAtomDerivatives(0,i) ); - } - setBoxDerivatives(valuex,cvout.virial[0]); - valuex->set(value[0]); - - for(unsigned i=0; i<2; ++i) { - setAtomsDerivatives(valuey,i,cvout.getAtomDerivatives(1,i) ); - } - setBoxDerivatives(valuey,cvout.virial[1]); - valuey->set(value[1]); - - for(unsigned i=0; i<2; ++i) { - setAtomsDerivatives(valuez,i,cvout.getAtomDerivatives(2,i) ); - } - setBoxDerivatives(valuez,cvout.virial[2]); - valuez->set(value[2]); - } else if( scaled_components ) { - ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) ); - calculateCV( ColvarInput::createColvarInput( 2, getPositions(), this ), cvout ); - - Value* valuea=getPntrToComponent("a"); - Value* valueb=getPntrToComponent("b"); - Value* valuec=getPntrToComponent("c"); - for(unsigned i=0; i<2; ++i) { - setAtomsDerivatives(valuea,i,cvout.getAtomDerivatives(0,i) ); - } - valuea->set(value[0]); - for(unsigned i=0; i<2; ++i) { - setAtomsDerivatives(valueb,i,cvout.getAtomDerivatives(1,i) ); - } - valueb->set(value[1]); - for(unsigned i=0; i<2; ++i) { - setAtomsDerivatives(valuec,i,cvout.getAtomDerivatives(2,i) ); - } - valuec->set(value[2]); - } else { - ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) ); - calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout ); - for(unsigned i=0; i<2; ++i) { - setAtomsDerivatives(i,cvout.getAtomDerivatives(0,i) ); - } - setBoxDerivatives(cvout.virial[0]); - setValue (value[0]); - } -} - -void Distance::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) { - Vector distance=delta(cvin.pos[0],cvin.pos[1]); - - if(cvin.mode==1) { - cvout.derivs[0][0] = Vector(-1,0,0); - cvout.derivs[0][1] = Vector(+1,0,0); - cvout.values[0] = distance[0]; - - cvout.derivs[1][0] = Vector(0,-1,0); - cvout.derivs[1][1] = Vector(0,+1,0); - cvout.values[1] = distance[1]; - - cvout.derivs[2][0] = Vector(0,0,-1); - cvout.derivs[2][1] = Vector(0,0,+1); - cvout.values[2] = distance[2]; - ColvarInput::setBoxDerivativesNoPbc( cvin, cvout ); - } else if(cvin.mode==2) { - Vector d=cvin.pbc.realToScaled(distance); - cvout.derivs[0][0] = matmul(cvin.pbc.getInvBox(),Vector(-1,0,0)); - cvout.derivs[0][1] = matmul(cvin.pbc.getInvBox(),Vector(+1,0,0)); - cvout.values[0] = Tools::pbc(d[0]); - cvout.derivs[1][0] = matmul(cvin.pbc.getInvBox(),Vector(0,-1,0)); - cvout.derivs[1][1] = matmul(cvin.pbc.getInvBox(),Vector(0,+1,0)); - cvout.values[1] = Tools::pbc(d[1]); - cvout.derivs[2][0] = matmul(cvin.pbc.getInvBox(),Vector(0,0,-1)); - cvout.derivs[2][1] = matmul(cvin.pbc.getInvBox(),Vector(0,0,+1)); - cvout.values[2] = Tools::pbc(d[2]); - } else { - const double value=distance.modulo(); - const double invvalue=1.0/value; - cvout.derivs[0][0] = -invvalue*distance; - cvout.derivs[0][1] = invvalue*distance; - ColvarInput::setBoxDerivativesNoPbc( cvin, cvout ); - cvout.values[0] = value; - } -} - -} -} - - - +} //namespace colvar +} //namespace PLMD diff --git a/src/colvar/Distance.h b/src/colvar/Distance.h new file mode 100644 index 0000000000..cd4eaaef87 --- /dev/null +++ b/src/colvar/Distance.h @@ -0,0 +1,249 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2011-2023 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see . ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#ifndef __PLUMED_colvar_Distance_h +#define __PLUMED_colvar_Distance_h +#include "Colvar.h" +#include "MultiColvarTemplate.h" +#include "tools/Pbc.h" + +namespace PLMD { +namespace colvar { + +class Distance : public Colvar { + bool components; + bool scaled_components; + bool pbc; + + std::vector value; + std::vector derivs; +public: + static void registerKeywords( Keywords& keys ); + explicit Distance(const ActionOptions&); + static void parseAtomList( const int& num, + std::vector& t, + ActionAtomistic* aa ); + static unsigned getModeAndSetupValues( ActionWithValue* av ); +// active methods: + void calculate() override; + static void calculateCV( const ColvarInput& cvin, + ColvarOutput& cvout ); +}; + +void Distance::registerKeywords( Keywords& keys ) { + Colvar::registerKeywords( keys ); + keys.setDisplayName("DISTANCE"); + constexpr auto scalarOrVector = Keywords::componentType::scalar + | Keywords::componentType::vector; + keys.add("atoms","ATOMS","the pair of atom that we are calculating the distance between"); + keys.addFlag("COMPONENTS",false,"calculate the x, y and z components of the distance separately and store them as label.x, label.y and label.z"); + keys.addFlag("SCALED_COMPONENTS",false,"calculate the a, b and c scaled components of the distance separately and store them as label.a, label.b and label.c"); + keys.addOutputComponent("x","COMPONENTS",scalarOrVector,"the x-component of the vector connecting the two atoms"); + keys.addOutputComponent("y","COMPONENTS",scalarOrVector,"the y-component of the vector connecting the two atoms"); + keys.addOutputComponent("z","COMPONENTS",scalarOrVector,"the z-component of the vector connecting the two atoms"); + keys.addOutputComponent("a","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the first lattice vector of the vector connecting the two atoms"); + keys.addOutputComponent("b","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the second lattice vector of the vector connecting the two atoms"); + keys.addOutputComponent("c","SCALED_COMPONENTS",scalarOrVector,"the normalized projection on the third lattice vector of the vector connecting the two atoms"); + keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log"); + keys.setValueDescription(scalarOrVector,"the DISTANCE between this pair of atoms"); + keys.addDOI("10.1007/978-1-4939-9608-7_21"); + keys.reset_style("NUMERICAL_DERIVATIVES","hidden"); +} + +Distance::Distance(const ActionOptions&ao): + PLUMED_COLVAR_INIT(ao), + components(false), + scaled_components(false), + pbc(true), + value(1) { + std::vector atoms; + parseAtomList(-1,atoms,this); + if(atoms.size()!=2) { + error("Number of specified atoms should be 2"); + } + + bool nopbc=!pbc; + parseFlag("NOPBC",nopbc); + pbc=!nopbc; + + if(pbc) { + log.printf(" using periodic boundary conditions\n"); + } else { + log.printf(" without periodic boundary conditions\n"); + } + + unsigned mode = getModeAndSetupValues( this ); + if(mode==1) { + components=true; + } else if(mode==2) { + scaled_components=true; + } + if( components || scaled_components ) { + value.resize(3); + } + requestAtoms(atoms); +} + +void Distance::parseAtomList( const int& num, + std::vector& t, + ActionAtomistic* aa ) { + aa->parseAtomList("ATOMS",num,t); + if( t.size()==2 ) { + aa->log.printf(" between atoms %d %d\n",t[0].serial(),t[1].serial()); + } +} + +unsigned Distance::getModeAndSetupValues( ActionWithValue* av ) { + bool c; + av->parseFlag("COMPONENTS",c); + bool sc; + av->parseFlag("SCALED_COMPONENTS",sc); + if( c && sc ) { + av->error("COMPONENTS and SCALED_COMPONENTS are not compatible"); + } + + if(c) { + av->addComponentWithDerivatives("x"); + av->componentIsNotPeriodic("x"); + av->addComponentWithDerivatives("y"); + av->componentIsNotPeriodic("y"); + av->addComponentWithDerivatives("z"); + av->componentIsNotPeriodic("z"); + av->log<<" WARNING: components will not have the proper periodicity - see manual\n"; + return 1; + } else if(sc) { + av->addComponentWithDerivatives("a"); + av->componentIsPeriodic("a","-0.5","+0.5"); + av->addComponentWithDerivatives("b"); + av->componentIsPeriodic("b","-0.5","+0.5"); + av->addComponentWithDerivatives("c"); + av->componentIsPeriodic("c","-0.5","+0.5"); + return 2; + } + av->addValueWithDerivatives(); + av->setNotPeriodic(); + return 0; +} + +// calculator +void Distance::calculate() { + + if(pbc) { + makeWhole(); + } + + if( components ) { + ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) ); + calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), cvout ); + Value* valuex=getPntrToComponent("x"); + Value* valuey=getPntrToComponent("y"); + Value* valuez=getPntrToComponent("z"); + + for(unsigned i=0; i<2; ++i) { + setAtomsDerivatives(valuex,i,cvout.getAtomDerivatives(0,i) ); + } + setBoxDerivatives(valuex,cvout.virial[0]); + valuex->set(value[0]); + + for(unsigned i=0; i<2; ++i) { + setAtomsDerivatives(valuey,i,cvout.getAtomDerivatives(1,i) ); + } + setBoxDerivatives(valuey,cvout.virial[1]); + valuey->set(value[1]); + + for(unsigned i=0; i<2; ++i) { + setAtomsDerivatives(valuez,i,cvout.getAtomDerivatives(2,i) ); + } + setBoxDerivatives(valuez,cvout.virial[2]); + valuez->set(value[2]); + } else if( scaled_components ) { + ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) ); + calculateCV( ColvarInput::createColvarInput( 2, getPositions(), this ), cvout ); + + Value* valuea=getPntrToComponent("a"); + Value* valueb=getPntrToComponent("b"); + Value* valuec=getPntrToComponent("c"); + for(unsigned i=0; i<2; ++i) { + setAtomsDerivatives(valuea,i,cvout.getAtomDerivatives(0,i) ); + } + valuea->set(value[0]); + for(unsigned i=0; i<2; ++i) { + setAtomsDerivatives(valueb,i,cvout.getAtomDerivatives(1,i) ); + } + valueb->set(value[1]); + for(unsigned i=0; i<2; ++i) { + setAtomsDerivatives(valuec,i,cvout.getAtomDerivatives(2,i) ); + } + valuec->set(value[2]); + } else { + ColvarOutput cvout( ColvarOutput::createColvarOutput(value, derivs, this) ); + calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout ); + for(unsigned i=0; i<2; ++i) { + setAtomsDerivatives(i,cvout.getAtomDerivatives(0,i) ); + } + setBoxDerivatives(cvout.virial[0]); + setValue (value[0]); + } +} + +void Distance::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) { + Vector distance=delta(cvin.pos[0],cvin.pos[1]); + + if(cvin.mode==1) { + cvout.derivs[0][0] = Vector(-1,0,0); + cvout.derivs[0][1] = Vector(+1,0,0); + cvout.values[0] = distance[0]; + + cvout.derivs[1][0] = Vector(0,-1,0); + cvout.derivs[1][1] = Vector(0,+1,0); + cvout.values[1] = distance[1]; + + cvout.derivs[2][0] = Vector(0,0,-1); + cvout.derivs[2][1] = Vector(0,0,+1); + cvout.values[2] = distance[2]; + ColvarInput::setBoxDerivativesNoPbc( cvin, cvout ); + } else if(cvin.mode==2) { + Vector d=cvin.pbc.realToScaled(distance); + cvout.derivs[0][0] = matmul(cvin.pbc.getInvBox(),Vector(-1,0,0)); + cvout.derivs[0][1] = matmul(cvin.pbc.getInvBox(),Vector(+1,0,0)); + cvout.values[0] = Tools::pbc(d[0]); + cvout.derivs[1][0] = matmul(cvin.pbc.getInvBox(),Vector(0,-1,0)); + cvout.derivs[1][1] = matmul(cvin.pbc.getInvBox(),Vector(0,+1,0)); + cvout.values[1] = Tools::pbc(d[1]); + cvout.derivs[2][0] = matmul(cvin.pbc.getInvBox(),Vector(0,0,-1)); + cvout.derivs[2][1] = matmul(cvin.pbc.getInvBox(),Vector(0,0,+1)); + cvout.values[2] = Tools::pbc(d[2]); + } else { + const double value=distance.modulo(); + const double invvalue=1.0/value; + cvout.derivs[0][0] = -invvalue*distance; + cvout.derivs[0][1] = invvalue*distance; + ColvarInput::setBoxDerivativesNoPbc( cvin, cvout ); + cvout.values[0] = value; + } +} + +} //namespace colvar +} //namespace PLMD + +#endif //__PLUMED_colvar_Distance_h + + diff --git a/src/colvar/MultiColvarTemplate.h b/src/colvar/MultiColvarTemplate.h index b8db6885fc..82e5ea4b7b 100644 --- a/src/colvar/MultiColvarTemplate.h +++ b/src/colvar/MultiColvarTemplate.h @@ -47,11 +47,12 @@ struct MultiColvarInput { } }; -template +template class MultiColvarTemplate : public ActionWithVector { public: using input_type = MultiColvarInput; - using PTM = ParallelTaskManager>; + using mytype=MultiColvarTemplate; + using PTM =typename myPTM::template PTM; constexpr static size_t virialSize = 9; private: /// The parallel task manager @@ -89,8 +90,8 @@ class MultiColvarTemplate : public ActionWithVector { ForceIndexHolder force_indices ); }; -template -void MultiColvarTemplate::registerKeywords(Keywords& keys ) { +template +void MultiColvarTemplate::registerKeywords(Keywords& keys ) { T::registerKeywords( keys ); PTM::registerKeywords( keys ); keys.addInputKeyword("optional","MASK","vector","the label for a sparse vector that should be used to determine which elements of the vector should be computed"); @@ -105,8 +106,8 @@ void MultiColvarTemplate::registerKeywords(Keywords& keys ) { } } -template -MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): +template +MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): Action(ao), ActionWithVector(ao), taskmanager(this), @@ -170,40 +171,40 @@ MultiColvarTemplate::MultiColvarTemplate(const ActionOptions&ao): taskmanager.setActionInput( MultiColvarInput{ usepbc, mode, natoms_per_task }); } -template -unsigned MultiColvarTemplate::getNumberOfDerivatives() { +template +unsigned MultiColvarTemplate::getNumberOfDerivatives() { return 3*getNumberOfAtoms()+9; } -template -void MultiColvarTemplate::calculate() { +template +void MultiColvarTemplate::calculate() { if( wholemolecules ) { makeWhole(); } taskmanager.runAllTasks(); } -template -void MultiColvarTemplate::applyNonZeroRankForces( std::vector& outforces ) { +template +void MultiColvarTemplate::applyNonZeroRankForces( std::vector& outforces ) { taskmanager.applyForces( outforces ); } -template -void MultiColvarTemplate::addValueWithDerivatives( const std::vector& shape ) { +template +void MultiColvarTemplate::addValueWithDerivatives( const std::vector& shape ) { std::vector s(1); s[0]=getNumberOfAtoms() / natoms_per_task; addValue( s ); } -template -void MultiColvarTemplate::addComponentWithDerivatives( const std::string& name, const std::vector& shape ) { +template +void MultiColvarTemplate::addComponentWithDerivatives( const std::string& name, const std::vector& shape ) { std::vector s(1); s[0]=getNumberOfAtoms() / natoms_per_task; addComponent( name, s ); } -template -void MultiColvarTemplate::getInputData( std::vector& inputdata ) const { +template +void MultiColvarTemplate::getInputData( std::vector& inputdata ) const { std::size_t ntasks = getConstPntrToComponent(0)->getNumberOfStoredValues(); if( inputdata.size()!=5*natoms_per_task*ntasks ) { inputdata.resize( 5*natoms_per_task*ntasks ); @@ -231,8 +232,8 @@ void MultiColvarTemplate::getInputData( std::vector& inputdata ) cons } } -template -void MultiColvarTemplate::performTask( std::size_t task_index, +template +void MultiColvarTemplate::performTask( std::size_t task_index, const MultiColvarInput& actiondata, ParallelActionsInput& input, ParallelActionsOutput& output ) { @@ -293,14 +294,14 @@ void MultiColvarTemplate::performTask( std::size_t task_index, cvout ); } -template -int MultiColvarTemplate::getNumberOfValuesPerTask( std::size_t task_index, +template +int MultiColvarTemplate::getNumberOfValuesPerTask( std::size_t task_index, const MultiColvarInput& actiondata ) { return 1; } -template -void MultiColvarTemplate::getForceIndices( std::size_t task_index, +template +void MultiColvarTemplate::getForceIndices( std::size_t task_index, std::size_t colno, std::size_t ntotal_force, const MultiColvarInput& actiondata, diff --git a/src/colvar/Plane.cpp b/src/colvar/Plane.cpp index a88f2ea845..1a4b48ecc5 100644 --- a/src/colvar/Plane.cpp +++ b/src/colvar/Plane.cpp @@ -19,18 +19,12 @@ You should have received a copy of the GNU Lesser General Public License along with plumed. If not, see . +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#ifdef __PLUMED_HAS_OPENACC -#define __PLUMED_USE_OPENACC 1 -#endif //__PLUMED_HAS_OPENACC -#include "Colvar.h" +#include "Plane.h" #include "ColvarShortcut.h" #include "core/ActionRegister.h" #include "MultiColvarTemplate.h" #include "tools/Pbc.h" -#include -#include - //+PLUMEDOC COLVAR PLANE /* Calculate the plane perpendicular to two vectors in order to represent the orientation of a planar molecule. @@ -78,126 +72,13 @@ nucleation of molecular crystals such as [SMAC](SMAC.md). namespace PLMD { namespace colvar { -class Plane : public Colvar { -private: - bool pbc; - std::vector value; - std::vector derivs; -public: - static void registerKeywords( Keywords& keys ); - explicit Plane(const ActionOptions&); - static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); -// active methods: - void calculate() override; - static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ); -}; - typedef ColvarShortcut PlaneShortcut; PLUMED_REGISTER_ACTION(PlaneShortcut,"PLANE") PLUMED_REGISTER_ACTION(Plane,"PLANE_SCALAR") typedef MultiColvarTemplate PlaneMulti; PLUMED_REGISTER_ACTION(PlaneMulti,"PLANE_VECTOR") -void Plane::registerKeywords( Keywords& keys ) { - Colvar::registerKeywords( keys ); - keys.setDisplayName("PLANE"); - keys.add("atoms","ATOMS","the three or four atoms whose plane we are computing"); - keys.addOutputComponent("x","default","scalar/vector","the x-component of the vector that is normal to the plane containing the atoms"); - keys.addOutputComponent("y","default","scalar/vector","the y-component of the vector that is normal to the plane containing the atoms"); - keys.addOutputComponent("z","default","scalar/vector","the z-component of the vector that is normal to the plane containing the atoms"); - keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log"); - keys.reset_style("NUMERICAL_DERIVATIVES","hidden"); -} - -void Plane::parseAtomList( const int& num, std::vector& atoms, ActionAtomistic* aa ) { - aa->parseAtomList("ATOMS",num,atoms); - if(atoms.size()==3) { - aa->log.printf(" containing atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial()); - atoms.resize(4); - atoms[3]=atoms[2]; - atoms[2]=atoms[1]; - } else if(atoms.size()==4) { - aa->log.printf(" containing lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial()); - } else if( num<0 || atoms.size()>0 ) { - aa->error("Number of specified atoms should be either 3 or 4"); - } -} - -unsigned Plane::getModeAndSetupValues( ActionWithValue* av ) { - av->addComponentWithDerivatives("x"); - av->componentIsNotPeriodic("x"); - av->addComponentWithDerivatives("y"); - av->componentIsNotPeriodic("y"); - av->addComponentWithDerivatives("z"); - av->componentIsNotPeriodic("z"); - return 0; -} -Plane::Plane(const ActionOptions&ao): - PLUMED_COLVAR_INIT(ao), - pbc(true), - value(3) { - std::vector atoms; - parseAtomList(-1,atoms,this); - bool nopbc=!pbc; - parseFlag("NOPBC",nopbc); - pbc=!nopbc; - - if(pbc) { - log.printf(" using periodic boundary conditions\n"); - } else { - log.printf(" without periodic boundary conditions\n"); - } - - unsigned mode = getModeAndSetupValues( this ); - requestAtoms(atoms); - checkRead(); } - -void Plane::calculate() { - - if(pbc) { - makeWhole(); - } - ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this); - calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout ); - setValue( value[0] ); - for(unsigned i=0; i. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#ifndef __PLUMED_colvar_Plane_h +#define __PLUMED_colvar_Plane_h +#include "Colvar.h" +#include "ColvarInput.h" +#include "tools/Pbc.h" + +#include +#include + +namespace PLMD { +namespace colvar { + +class Plane : public Colvar { +private: + bool pbc; + std::vector value; + std::vector derivs; +public: + static void registerKeywords( Keywords& keys ); + explicit Plane(const ActionOptions&); + static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); + static unsigned getModeAndSetupValues( ActionWithValue* av ); +// active methods: + void calculate() override; + static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ); +}; + +void Plane::registerKeywords( Keywords& keys ) { + Colvar::registerKeywords( keys ); + keys.setDisplayName("PLANE"); + keys.add("atoms","ATOMS","the three or four atoms whose plane we are computing"); + keys.addOutputComponent("x","default","scalar/vector","the x-component of the vector that is normal to the plane containing the atoms"); + keys.addOutputComponent("y","default","scalar/vector","the y-component of the vector that is normal to the plane containing the atoms"); + keys.addOutputComponent("z","default","scalar/vector","the z-component of the vector that is normal to the plane containing the atoms"); + keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log"); + keys.reset_style("NUMERICAL_DERIVATIVES","hidden"); +} + +void Plane::parseAtomList( const int& num, std::vector& atoms, ActionAtomistic* aa ) { + aa->parseAtomList("ATOMS",num,atoms); + if(atoms.size()==3) { + aa->log.printf(" containing atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial()); + atoms.resize(4); + atoms[3]=atoms[2]; + atoms[2]=atoms[1]; + } else if(atoms.size()==4) { + aa->log.printf(" containing lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial()); + } else if( num<0 || atoms.size()>0 ) { + aa->error("Number of specified atoms should be either 3 or 4"); + } +} + +unsigned Plane::getModeAndSetupValues( ActionWithValue* av ) { + av->addComponentWithDerivatives("x"); + av->componentIsNotPeriodic("x"); + av->addComponentWithDerivatives("y"); + av->componentIsNotPeriodic("y"); + av->addComponentWithDerivatives("z"); + av->componentIsNotPeriodic("z"); + return 0; +} + +Plane::Plane(const ActionOptions&ao): + PLUMED_COLVAR_INIT(ao), + pbc(true), + value(3) { + std::vector atoms; + parseAtomList(-1,atoms,this); + bool nopbc=!pbc; + parseFlag("NOPBC",nopbc); + pbc=!nopbc; + + if(pbc) { + log.printf(" using periodic boundary conditions\n"); + } else { + log.printf(" without periodic boundary conditions\n"); + } + + unsigned mode = getModeAndSetupValues( this ); + requestAtoms(atoms); + checkRead(); +} + +void Plane::calculate() { + + if(pbc) { + makeWhole(); + } + ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this); + calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout ); + setValue( value[0] ); + for(unsigned i=0; i. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#ifdef __PLUMED_HAS_OPENACC -#define __PLUMED_USE_OPENACC 1 -#endif //__PLUMED_HAS_OPENACC -#include "Colvar.h" +#include "Position.h" #include "ColvarShortcut.h" #include "MultiColvarTemplate.h" #include "core/ActionRegister.h" -#include "tools/Pbc.h" namespace PLMD { namespace colvar { @@ -106,164 +102,12 @@ PRINT ARG=p.x,p.y,p.z */ //+ENDPLUMEDOC -class Position : public Colvar { - bool scaled_components; - bool pbc; - std::vector value; - std::vector derivs; -public: - static void registerKeywords( Keywords& keys ); - explicit Position(const ActionOptions&); - static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); -// active methods: - void calculate() override; - static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ); -}; - typedef ColvarShortcut PositionShortcut; PLUMED_REGISTER_ACTION(PositionShortcut,"POSITION") PLUMED_REGISTER_ACTION(Position,"POSITION_SCALAR") typedef MultiColvarTemplate PositionMulti; PLUMED_REGISTER_ACTION(PositionMulti,"POSITION_VECTOR") -void Position::registerKeywords( Keywords& keys ) { - Colvar::registerKeywords( keys ); - keys.setDisplayName("POSITION"); - keys.add("atoms","ATOM","the atom number"); - keys.add("atoms","ATOMS","the atom numbers that you would like to use the positions of"); - keys.addFlag("WHOLEMOLECULES",false,"if this is a vector of positions do you want to make the positions into a whole before"); - keys.addFlag("SCALED_COMPONENTS",false,"calculate the a, b and c scaled components of the position separately and store them as label.a, label.b and label.c"); - keys.addOutputComponent("x","default","scalar/vector","the x-component of the atom position"); - keys.addOutputComponent("y","default","scalar/vector","the y-component of the atom position"); - keys.addOutputComponent("z","default","scalar/vector","the z-component of the atom position"); - keys.addOutputComponent("a","SCALED_COMPONENTS","scalar/vector","the normalized projection on the first lattice vector of the atom position"); - keys.addOutputComponent("b","SCALED_COMPONENTS","scalar/vector","the normalized projection on the second lattice vector of the atom position"); - keys.addOutputComponent("c","SCALED_COMPONENTS","scalar/vector","the normalized projection on the third lattice vector of the atom position"); - keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log"); - keys.reset_style("NUMERICAL_DERIVATIVES","hidden"); -} - -Position::Position(const ActionOptions&ao): - PLUMED_COLVAR_INIT(ao), - scaled_components(false), - pbc(true), - value(3) { - std::vector atoms; - parseAtomList(-1,atoms,this); - unsigned mode=getModeAndSetupValues(this); - if( mode==1 ) { - scaled_components=true; - } - - bool nopbc=!pbc; - parseFlag("NOPBC",nopbc); - pbc=!nopbc; - checkRead(); - - if(pbc) { - log.printf(" using periodic boundary conditions\n"); - } else { - log.printf(" without periodic boundary conditions\n"); - } - - requestAtoms(atoms); -} - -void Position::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { - aa->parseAtomList("ATOM",num,t); - if( t.size()==1 ) { - aa->log.printf(" for atom %d\n",t[0].serial()); - } else if( num<0 || t.size()!=0 ) { - aa->error("Number of specified atoms should be 1"); - } -} - -unsigned Position::getModeAndSetupValues( ActionWithValue* av ) { - bool sc; - av->parseFlag("SCALED_COMPONENTS",sc); - if(sc) { - av->addComponentWithDerivatives("a"); - av->componentIsPeriodic("a","-0.5","+0.5"); - av->addComponentWithDerivatives("b"); - av->componentIsPeriodic("b","-0.5","+0.5"); - av->addComponentWithDerivatives("c"); - av->componentIsPeriodic("c","-0.5","+0.5"); - return 1; - } - av->addComponentWithDerivatives("x"); - av->componentIsNotPeriodic("x"); - av->addComponentWithDerivatives("y"); - av->componentIsNotPeriodic("y"); - av->addComponentWithDerivatives("z"); - av->componentIsNotPeriodic("z"); - av->log<<" WARNING: components will not have the proper periodicity - see manual\n"; - return 0; -} - -// calculator -void Position::calculate() { - - std::vector distance(1); - if(pbc) { - distance[0]=pbcDistance(Vector(0.0,0.0,0.0),getPosition(0)); - } else { - distance[0]=delta(Vector(0.0,0.0,0.0),getPosition(0)); - } - - ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this); - if(scaled_components) { - calculateCV( ColvarInput::createColvarInput( 1, distance, this ), cvout ); - Value* valuea=getPntrToComponent("a"); - Value* valueb=getPntrToComponent("b"); - Value* valuec=getPntrToComponent("c"); - setAtomsDerivatives (valuea,0,cvout.getAtomDerivatives(0,0)); - valuea->set(value[0]); - setAtomsDerivatives (valueb,0,cvout.getAtomDerivatives(1,0)); - valueb->set(value[1]); - setAtomsDerivatives (valuec,0,cvout.getAtomDerivatives(2,0)); - valuec->set(value[2]); - } else { - calculateCV( ColvarInput::createColvarInput( 0, distance, this ), cvout ); - Value* valuex=getPntrToComponent("x"); - Value* valuey=getPntrToComponent("y"); - Value* valuez=getPntrToComponent("z"); - - setAtomsDerivatives (valuex,0,cvout.getAtomDerivatives(0,0)); - setBoxDerivatives (valuex,cvout.virial[0]); - valuex->set(value[0]); - - setAtomsDerivatives (valuey,0,cvout.getAtomDerivatives(1,0)); - setBoxDerivatives (valuey,cvout.virial[1]); - valuey->set(value[1]); - - setAtomsDerivatives (valuez,0,cvout.getAtomDerivatives(2,0)); - setBoxDerivatives (valuez,cvout.virial[2]); - valuez->set(value[2]); - } -} - -void Position::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) { - if( cvin.mode==1 ) { - Vector d=cvin.pbc.realToScaled(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2])); - cvout.values[0]=Tools::pbc(d[0]); - cvout.values[1]=Tools::pbc(d[1]); - cvout.values[2]=Tools::pbc(d[2]); - cvout.derivs[0][0]=matmul(cvin.pbc.getInvBox(),Vector(+1,0,0)); - cvout.derivs[1][0]=matmul(cvin.pbc.getInvBox(),Vector(0,+1,0)); - cvout.derivs[2][0]=matmul(cvin.pbc.getInvBox(),Vector(0,0,+1)); - } else { - for(unsigned i=0; i<3; ++i) { - cvout.values[i]=cvin.pos[0][i]; - } - cvout.derivs[0][0]=Vector(+1,0,0); - cvout.derivs[1][0]=Vector(0,+1,0); - cvout.derivs[2][0]=Vector(0,0,+1); - cvout.virial.set(0, Tensor(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]),Vector(-1,0,0)) ); - cvout.virial.set(1, Tensor(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]),Vector(0,-1,0)) ); - cvout.virial.set(2, Tensor(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]),Vector(0,0,-1)) ); - } -} } } diff --git a/src/colvar/Position.h b/src/colvar/Position.h new file mode 100644 index 0000000000..ce56cdc7a8 --- /dev/null +++ b/src/colvar/Position.h @@ -0,0 +1,186 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2014-2023 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see . ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#ifndef __PLUMED_colvar_Position_h +#define __PLUMED_colvar_Position_h +#include "Colvar.h" +#include "tools/Pbc.h" +#include "ColvarInput.h" + +namespace PLMD { +namespace colvar { + +class Position : public Colvar { + bool scaled_components; + bool pbc; + std::vector value; + std::vector derivs; +public: + static void registerKeywords( Keywords& keys ); + explicit Position(const ActionOptions&); + static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); + static unsigned getModeAndSetupValues( ActionWithValue* av ); +// active methods: + void calculate() override; + static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ); +}; + +void Position::registerKeywords( Keywords& keys ) { + Colvar::registerKeywords( keys ); + keys.setDisplayName("POSITION"); + keys.add("atoms","ATOM","the atom number"); + keys.add("atoms","ATOMS","the atom numbers that you would like to use the positions of"); + keys.addFlag("WHOLEMOLECULES",false,"if this is a vector of positions do you want to make the positions into a whole before"); + keys.addFlag("SCALED_COMPONENTS",false,"calculate the a, b and c scaled components of the position separately and store them as label.a, label.b and label.c"); + keys.addOutputComponent("x","default","scalar/vector","the x-component of the atom position"); + keys.addOutputComponent("y","default","scalar/vector","the y-component of the atom position"); + keys.addOutputComponent("z","default","scalar/vector","the z-component of the atom position"); + keys.addOutputComponent("a","SCALED_COMPONENTS","scalar/vector","the normalized projection on the first lattice vector of the atom position"); + keys.addOutputComponent("b","SCALED_COMPONENTS","scalar/vector","the normalized projection on the second lattice vector of the atom position"); + keys.addOutputComponent("c","SCALED_COMPONENTS","scalar/vector","the normalized projection on the third lattice vector of the atom position"); + keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log"); + keys.reset_style("NUMERICAL_DERIVATIVES","hidden"); +} + +Position::Position(const ActionOptions&ao): + PLUMED_COLVAR_INIT(ao), + scaled_components(false), + pbc(true), + value(3) { + std::vector atoms; + parseAtomList(-1,atoms,this); + unsigned mode=getModeAndSetupValues(this); + if( mode==1 ) { + scaled_components=true; + } + + bool nopbc=!pbc; + parseFlag("NOPBC",nopbc); + pbc=!nopbc; + checkRead(); + + if(pbc) { + log.printf(" using periodic boundary conditions\n"); + } else { + log.printf(" without periodic boundary conditions\n"); + } + + requestAtoms(atoms); +} + +void Position::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { + aa->parseAtomList("ATOM",num,t); + if( t.size()==1 ) { + aa->log.printf(" for atom %d\n",t[0].serial()); + } else if( num<0 || t.size()!=0 ) { + aa->error("Number of specified atoms should be 1"); + } +} + +unsigned Position::getModeAndSetupValues( ActionWithValue* av ) { + bool sc; + av->parseFlag("SCALED_COMPONENTS",sc); + if(sc) { + av->addComponentWithDerivatives("a"); + av->componentIsPeriodic("a","-0.5","+0.5"); + av->addComponentWithDerivatives("b"); + av->componentIsPeriodic("b","-0.5","+0.5"); + av->addComponentWithDerivatives("c"); + av->componentIsPeriodic("c","-0.5","+0.5"); + return 1; + } + av->addComponentWithDerivatives("x"); + av->componentIsNotPeriodic("x"); + av->addComponentWithDerivatives("y"); + av->componentIsNotPeriodic("y"); + av->addComponentWithDerivatives("z"); + av->componentIsNotPeriodic("z"); + av->log<<" WARNING: components will not have the proper periodicity - see manual\n"; + return 0; +} + +// calculator +void Position::calculate() { + + std::vector distance(1); + if(pbc) { + distance[0]=pbcDistance(Vector(0.0,0.0,0.0),getPosition(0)); + } else { + distance[0]=delta(Vector(0.0,0.0,0.0),getPosition(0)); + } + + ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this); + if(scaled_components) { + calculateCV( ColvarInput::createColvarInput( 1, distance, this ), cvout ); + Value* valuea=getPntrToComponent("a"); + Value* valueb=getPntrToComponent("b"); + Value* valuec=getPntrToComponent("c"); + setAtomsDerivatives (valuea,0,cvout.getAtomDerivatives(0,0)); + valuea->set(value[0]); + setAtomsDerivatives (valueb,0,cvout.getAtomDerivatives(1,0)); + valueb->set(value[1]); + setAtomsDerivatives (valuec,0,cvout.getAtomDerivatives(2,0)); + valuec->set(value[2]); + } else { + calculateCV( ColvarInput::createColvarInput( 0, distance, this ), cvout ); + Value* valuex=getPntrToComponent("x"); + Value* valuey=getPntrToComponent("y"); + Value* valuez=getPntrToComponent("z"); + + setAtomsDerivatives (valuex,0,cvout.getAtomDerivatives(0,0)); + setBoxDerivatives (valuex,cvout.virial[0]); + valuex->set(value[0]); + + setAtomsDerivatives (valuey,0,cvout.getAtomDerivatives(1,0)); + setBoxDerivatives (valuey,cvout.virial[1]); + valuey->set(value[1]); + + setAtomsDerivatives (valuez,0,cvout.getAtomDerivatives(2,0)); + setBoxDerivatives (valuez,cvout.virial[2]); + valuez->set(value[2]); + } +} + +void Position::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) { + if( cvin.mode==1 ) { + Vector d=cvin.pbc.realToScaled(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2])); + cvout.values[0]=Tools::pbc(d[0]); + cvout.values[1]=Tools::pbc(d[1]); + cvout.values[2]=Tools::pbc(d[2]); + cvout.derivs[0][0]=matmul(cvin.pbc.getInvBox(),Vector(+1,0,0)); + cvout.derivs[1][0]=matmul(cvin.pbc.getInvBox(),Vector(0,+1,0)); + cvout.derivs[2][0]=matmul(cvin.pbc.getInvBox(),Vector(0,0,+1)); + } else { + for(unsigned i=0; i<3; ++i) { + cvout.values[i]=cvin.pos[0][i]; + } + cvout.derivs[0][0]=Vector(+1,0,0); + cvout.derivs[1][0]=Vector(0,+1,0); + cvout.derivs[2][0]=Vector(0,0,+1); + cvout.virial.set(0, Tensor(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]),Vector(-1,0,0)) ); + cvout.virial.set(1, Tensor(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]),Vector(0,-1,0)) ); + cvout.virial.set(2, Tensor(Vector(cvin.pos[0][0],cvin.pos[0][1],cvin.pos[0][2]),Vector(0,0,-1)) ); + } +} + +} +} +#endif //__PLUMED_colvar_Position_h diff --git a/src/colvar/Torsion.cpp b/src/colvar/Torsion.cpp index a55cb3dedb..ca7cbddba5 100644 --- a/src/colvar/Torsion.cpp +++ b/src/colvar/Torsion.cpp @@ -19,15 +19,10 @@ You should have received a copy of the GNU Lesser General Public License along with plumed. If not, see . +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -//this is temporary: -#ifdef __PLUMED_HAS_OPENACC -#define __PLUMED_USE_OPENACC 1 -#endif //__PLUMED_HAS_OPENACC -#include "Colvar.h" +#include "Torsion.h" #include "ColvarShortcut.h" #include "MultiColvarTemplate.h" #include "core/ActionRegister.h" -#include "tools/Torsion.h" namespace PLMD { namespace colvar { @@ -115,185 +110,12 @@ PRINT ARG=t FILE=COLVAR */ //+ENDPLUMEDOC -class Torsion : public Colvar { - bool pbc; - bool do_cosine; - - std::vector value; - std::vector derivs; -public: - explicit Torsion(const ActionOptions&); - static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); - static unsigned getModeAndSetupValues( ActionWithValue* av ); -// active methods: - void calculate() override; - static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ); - static void registerKeywords(Keywords& keys); -}; - typedef ColvarShortcut TorsionShortcut; PLUMED_REGISTER_ACTION(TorsionShortcut,"TORSION") PLUMED_REGISTER_ACTION(Torsion,"TORSION_SCALAR") typedef MultiColvarTemplate TorsionMulti; PLUMED_REGISTER_ACTION(TorsionMulti,"TORSION_VECTOR") -void Torsion::registerKeywords(Keywords& keys) { - Colvar::registerKeywords( keys ); - keys.setDisplayName("TORSION"); - keys.add("atoms-1","ATOMS","the four atoms involved in the torsional angle"); - keys.add("atoms-2","AXIS","two atoms that define an axis. You can use this to find the angle in the plane perpendicular to the axis between the vectors specified using the VECTORA and VECTORB keywords."); - keys.add("atoms-2","VECTORA","two atoms that define a vector. You can use this in combination with VECTOR2 and AXIS"); - keys.add("atoms-2","VECTORB","two atoms that define a vector. You can use this in combination with VECTOR1 and AXIS"); - keys.addDeprecatedKeyword("VECTOR1","VECTORA"); - keys.addDeprecatedKeyword("VECTOR2","VECTORB"); - keys.addFlag("COSINE",false,"calculate cosine instead of dihedral"); - keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log"); - keys.setValueDescription("scalar/vector","the TORSION involving these atoms"); - keys.reset_style("NUMERICAL_DERIVATIVES","hidden"); -} - -Torsion::Torsion(const ActionOptions&ao): - PLUMED_COLVAR_INIT(ao), - pbc(true), - do_cosine(false), - value(1) { - std::vector atoms; - std::vector v1; - ActionAtomistic::parseAtomList("VECTOR1",v1); - if( v1.size()>0 ) { - std::vector v2; - ActionAtomistic::parseAtomList("VECTOR2",v2); - std::vector axis; - ActionAtomistic::parseAtomList("AXIS",axis); - if( !(v1.size()==2 && v2.size()==2 && axis.size()==2)) { - error("VECTOR1, VECTOR2 and AXIS should specify 2 atoms each"); - } - atoms.resize(6); - atoms[0]=v1[1]; - atoms[1]=v1[0]; - atoms[2]=axis[0]; - atoms[3]=axis[1]; - atoms[4]=v2[0]; - atoms[5]=v2[1]; - log.printf(" between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n", - v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial()); - } else { - parseAtomList(-1,atoms,this); - } - unsigned mode=getModeAndSetupValues(this); - if( mode==1 ) { - do_cosine=true; - } - - bool nopbc=!pbc; - parseFlag("NOPBC",nopbc); - pbc=!nopbc; - checkRead(); - - if(pbc) { - log.printf(" using periodic boundary conditions\n"); - } else { - log.printf(" without periodic boundary conditions\n"); - } - requestAtoms(atoms); -} - -void Torsion::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { - std::vector v1,v2,axis; - aa->parseAtomList("ATOMS",num,t); - aa->parseAtomList("VECTORA",num,v1); - aa->parseAtomList("VECTORB",num,v2); - aa->parseAtomList("AXIS",num,axis); - - if(t.size()==4) { - if(!(v1.empty() && v2.empty() && axis.empty())) { - aa->error("ATOMS keyword is not compatible with VECTORA, VECTORB and AXIS keywords"); - } - aa->log.printf(" between atoms %d %d %d %d\n",t[0].serial(),t[1].serial(),t[2].serial(),t[3].serial()); - t.resize(6); - t[5]=t[3]; - t[4]=t[2]; - t[3]=t[2]; - t[2]=t[1]; - } else if(t.empty()) { - if( num>0 && v1.empty() && v2.empty() && axis.empty() ) { - return; - } - if(!(v1.size()==2 && v2.size()==2 && axis.size()==2)) { - aa->error("VECTORA, VECTORB and AXIS should specify 2 atoms each"); - } - aa->log.printf(" between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n", - v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial()); - t.resize(6); - t[0]=v1[1]; - t[1]=v1[0]; - t[2]=axis[0]; - t[3]=axis[1]; - t[4]=v2[0]; - t[5]=v2[1]; - } else if( t.size()!=4 ) { - aa->error("ATOMS should specify 4 atoms"); - } -} - -unsigned Torsion::getModeAndSetupValues( ActionWithValue* av ) { - bool do_cos; - av->parseFlag("COSINE",do_cos); - if(do_cos) { - av->log.printf(" calculating cosine instead of torsion\n"); - } - - av->addValueWithDerivatives(); - if(!do_cos) { - av->setPeriodic("-pi","pi"); - return 0; - } - av->setNotPeriodic(); - return 1; -} - -// calculator -void Torsion::calculate() { - if(pbc) { - makeWhole(); - } - ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this); - if(do_cosine) { - calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), cvout ); - } else { - calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout ); - } - for(unsigned i=0; i<6; ++i) { - setAtomsDerivatives(i,cvout.getAtomDerivatives(0,i) ); - } - setValue(value[0]); - setBoxDerivatives( cvout.virial[0] ); -} - -void Torsion::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) { - Vector d0=delta(cvin.pos[1],cvin.pos[0]); - Vector d1=delta(cvin.pos[3],cvin.pos[2]); - Vector d2=delta(cvin.pos[5],cvin.pos[4]); - Vector dd0,dd1,dd2; - PLMD::Torsion t; - cvout.values[0] = t.compute(d0,d1,d2,dd0,dd1,dd2); - if(cvin.mode==1) { - dd0 *= -std::sin(cvout.values[0]); - dd1 *= -std::sin(cvout.values[0]); - dd2 *= -std::sin(cvout.values[0]); - cvout.values[0] = std::cos(cvout.values[0]); - } - cvout.derivs[0][0] = dd0; - cvout.derivs[0][1] = -dd0; - cvout.derivs[0][2] = dd1; - cvout.derivs[0][3] = -dd1; - cvout.derivs[0][4] = dd2; - cvout.derivs[0][5] = -dd2; - ColvarInput::setBoxDerivativesNoPbc( cvin, cvout ); -} - } } - - diff --git a/src/colvar/Torsion.h b/src/colvar/Torsion.h new file mode 100644 index 0000000000..00a2b9778d --- /dev/null +++ b/src/colvar/Torsion.h @@ -0,0 +1,204 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2011-2023 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see . ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#ifndef __PLUMED_colvar_Torsion_h +#define __PLUMED_colvar_Torsion_h +#include "Colvar.h" +#include "tools/Torsion.h" +#include "ColvarInput.h" + +namespace PLMD { +namespace colvar { + +class Torsion : public Colvar { + bool pbc; + bool do_cosine; + + std::vector value; + std::vector derivs; +public: + explicit Torsion(const ActionOptions&); + static void parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ); + static unsigned getModeAndSetupValues( ActionWithValue* av ); +// active methods: + void calculate() override; + static void calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ); + static void registerKeywords(Keywords& keys); +}; + +void Torsion::registerKeywords(Keywords& keys) { + Colvar::registerKeywords( keys ); + keys.setDisplayName("TORSION"); + keys.add("atoms-1","ATOMS","the four atoms involved in the torsional angle"); + keys.add("atoms-2","AXIS","two atoms that define an axis. You can use this to find the angle in the plane perpendicular to the axis between the vectors specified using the VECTORA and VECTORB keywords."); + keys.add("atoms-2","VECTORA","two atoms that define a vector. You can use this in combination with VECTOR2 and AXIS"); + keys.add("atoms-2","VECTORB","two atoms that define a vector. You can use this in combination with VECTOR1 and AXIS"); + keys.addDeprecatedKeyword("VECTOR1","VECTORA"); + keys.addDeprecatedKeyword("VECTOR2","VECTORB"); + keys.addFlag("COSINE",false,"calculate cosine instead of dihedral"); + keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log"); + keys.setValueDescription("scalar/vector","the TORSION involving these atoms"); + keys.reset_style("NUMERICAL_DERIVATIVES","hidden"); +} + +Torsion::Torsion(const ActionOptions&ao): + PLUMED_COLVAR_INIT(ao), + pbc(true), + do_cosine(false), + value(1) { + std::vector atoms; + std::vector v1; + ActionAtomistic::parseAtomList("VECTOR1",v1); + if( v1.size()>0 ) { + std::vector v2; + ActionAtomistic::parseAtomList("VECTOR2",v2); + std::vector axis; + ActionAtomistic::parseAtomList("AXIS",axis); + if( !(v1.size()==2 && v2.size()==2 && axis.size()==2)) { + error("VECTOR1, VECTOR2 and AXIS should specify 2 atoms each"); + } + atoms.resize(6); + atoms[0]=v1[1]; + atoms[1]=v1[0]; + atoms[2]=axis[0]; + atoms[3]=axis[1]; + atoms[4]=v2[0]; + atoms[5]=v2[1]; + log.printf(" between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n", + v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial()); + } else { + parseAtomList(-1,atoms,this); + } + unsigned mode=getModeAndSetupValues(this); + if( mode==1 ) { + do_cosine=true; + } + + bool nopbc=!pbc; + parseFlag("NOPBC",nopbc); + pbc=!nopbc; + checkRead(); + + if(pbc) { + log.printf(" using periodic boundary conditions\n"); + } else { + log.printf(" without periodic boundary conditions\n"); + } + requestAtoms(atoms); +} + +void Torsion::parseAtomList( const int& num, std::vector& t, ActionAtomistic* aa ) { + std::vector v1,v2,axis; + aa->parseAtomList("ATOMS",num,t); + aa->parseAtomList("VECTORA",num,v1); + aa->parseAtomList("VECTORB",num,v2); + aa->parseAtomList("AXIS",num,axis); + + if(t.size()==4) { + if(!(v1.empty() && v2.empty() && axis.empty())) { + aa->error("ATOMS keyword is not compatible with VECTORA, VECTORB and AXIS keywords"); + } + aa->log.printf(" between atoms %d %d %d %d\n",t[0].serial(),t[1].serial(),t[2].serial(),t[3].serial()); + t.resize(6); + t[5]=t[3]; + t[4]=t[2]; + t[3]=t[2]; + t[2]=t[1]; + } else if(t.empty()) { + if( num>0 && v1.empty() && v2.empty() && axis.empty() ) { + return; + } + if(!(v1.size()==2 && v2.size()==2 && axis.size()==2)) { + aa->error("VECTORA, VECTORB and AXIS should specify 2 atoms each"); + } + aa->log.printf(" between lines %d-%d and %d-%d, projected on the plane orthogonal to line %d-%d\n", + v1[0].serial(),v1[1].serial(),v2[0].serial(),v2[1].serial(),axis[0].serial(),axis[1].serial()); + t.resize(6); + t[0]=v1[1]; + t[1]=v1[0]; + t[2]=axis[0]; + t[3]=axis[1]; + t[4]=v2[0]; + t[5]=v2[1]; + } else if( t.size()!=4 ) { + aa->error("ATOMS should specify 4 atoms"); + } +} + +unsigned Torsion::getModeAndSetupValues( ActionWithValue* av ) { + bool do_cos; + av->parseFlag("COSINE",do_cos); + if(do_cos) { + av->log.printf(" calculating cosine instead of torsion\n"); + } + + av->addValueWithDerivatives(); + if(!do_cos) { + av->setPeriodic("-pi","pi"); + return 0; + } + av->setNotPeriodic(); + return 1; +} + +// calculator +void Torsion::calculate() { + if(pbc) { + makeWhole(); + } + ColvarOutput cvout = ColvarOutput::createColvarOutput(value,derivs,this); + if(do_cosine) { + calculateCV( ColvarInput::createColvarInput( 1, getPositions(), this ), cvout ); + } else { + calculateCV( ColvarInput::createColvarInput( 0, getPositions(), this ), cvout ); + } + for(unsigned i=0; i<6; ++i) { + setAtomsDerivatives(i,cvout.getAtomDerivatives(0,i) ); + } + setValue(value[0]); + setBoxDerivatives( cvout.virial[0] ); +} + +void Torsion::calculateCV( const ColvarInput& cvin, ColvarOutput& cvout ) { + Vector d0=delta(cvin.pos[1],cvin.pos[0]); + Vector d1=delta(cvin.pos[3],cvin.pos[2]); + Vector d2=delta(cvin.pos[5],cvin.pos[4]); + Vector dd0,dd1,dd2; + PLMD::Torsion t; + cvout.values[0] = t.compute(d0,d1,d2,dd0,dd1,dd2); + if(cvin.mode==1) { + dd0 *= -std::sin(cvout.values[0]); + dd1 *= -std::sin(cvout.values[0]); + dd2 *= -std::sin(cvout.values[0]); + cvout.values[0] = std::cos(cvout.values[0]); + } + cvout.derivs[0][0] = dd0; + cvout.derivs[0][1] = -dd0; + cvout.derivs[0][2] = dd1; + cvout.derivs[0][3] = -dd1; + cvout.derivs[0][4] = dd2; + cvout.derivs[0][5] = -dd2; + ColvarInput::setBoxDerivativesNoPbc( cvin, cvout ); +} + +} +} +#endif // __PLUMED_colvar_Torsion.h diff --git a/src/core/ActionShortcut.cpp b/src/core/ActionShortcut.cpp index ab03e9d01e..b694cd509a 100644 --- a/src/core/ActionShortcut.cpp +++ b/src/core/ActionShortcut.cpp @@ -149,7 +149,12 @@ void ActionShortcut::readInputLine( const std::string& input, bool saveline ) { Keywords thiskeys; actionRegister().getKeywords( av->getName(), thiskeys ); if( thiskeys.getDisplayName()!=getName() ) { - error("mismatch between display name of hidden action " + thiskeys.getDisplayName() + " and shortcut that creates it " + getName() ); + //NOTE: this was an error + // I deactivated it because it clashes with the gpu-openacc-by-LOAD does not load + // at moment of writing this I do not know if there is a way to to the check by passing the + // dlloader object to this function. + // bypassing this check makes my contraption work + warning("mismatch between display name of hidden action " + thiskeys.getDisplayName() + " and shortcut that creates it " + getName() ); } } } else { diff --git a/src/core/ParallelTaskManager.h b/src/core/ParallelTaskManager.h index b40bfaf5c9..d203e44acc 100644 --- a/src/core/ParallelTaskManager.h +++ b/src/core/ParallelTaskManager.h @@ -381,7 +381,7 @@ class ParallelTaskManager { // static constexpr bool has_custom_gather=PTMUtils::has_gatherForces_custom; // static constexpr bool has_GPU_gather=PTMUtils::has_gatherForces_GPU; // static constexpr size_t virialSize = PTMUtils::virialSize; -private: +protected: /// The underlying action for which we are managing parallel tasks ActionWithVector* action; /// The MPI communicator @@ -442,6 +442,11 @@ class ParallelTaskManager { View forces ); }; +struct defaultPTM { + template + using PTM=ParallelTaskManager; +}; + template void ParallelTaskManager::registerKeywords( Keywords& keys ) { keys.addFlag("USEGPU",false,"run this calculation on the GPU");