ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Update Colvars to version 2017-07-15 and support automated builds for it

Browse Source
This commit is contained in:
Giacomo Fiorin
2017-07-19 14:10:43 -04:00
parent cc9b3864bf
commit ee6cac826e
52 changed files with 3473 additions and 2413 deletions
Download Patch File Download Diff File Expand all files Collapse all files

BIN
doc/src/PDF/colvars-refman-lammps.pdf
View File

Binary file not shown.

102
lib/colvars/Install.py
View File

@ -1,27 +1,34 @@
#!/usr/bin/env python
# install.py tool to do a generic build of a library
# soft linked to by many of the lib/Install.py files
# used to automate the steps described in the corresponding lib/README
# Install.py tool to do automate build of Colvars
import sys,commands,os
from __future__ import print_function
import sys,os,subprocess
# help message
help = """
Syntax: python Install.py -m machine -e suffix
specify -m and optionally -e, order does not matter
Syntax from src dir: make lib-colvars args="-m machine -e suffix"
Syntax from lib/colvars dir: python Install.py -m machine -e suffix
specify -m and optionally -e, order does not matter
-m = peform a clean followed by "make -f Makefile.machine"
machine = suffix of a lib/Makefile.* file
machine = suffix of a lib/colvars/Makefile.* or of a
src/MAKE/MACHINES/Makefile.* file
-e = set EXTRAMAKE variable in Makefile.machine to Makefile.lammps.suffix
does not alter existing Makefile.machine
Examples:
make lib-colvars args="-m g++" # build COLVARS lib with GNU g++ compiler
"""
# print error message or help
def error(str=None):
if not str: print help
else: print "ERROR",str
if not str: print(help)
else: print("ERROR"),str
sys.exit()
# parse args
@ -31,17 +38,17 @@ nargs = len(args)
if nargs == 0: error()
machine = None
extraflag = 0
extraflag = False
iarg = 0
while iarg < nargs:
if args[iarg] == "-m":
if iarg+2 > nargs: error()
if iarg+2 > len(args): error()
machine = args[iarg+1]
iarg += 2
elif args[iarg] == "-e":
if iarg+2 > nargs: error()
extraflag = 1
if iarg+2 > len(args): error()
extraflag = True
suffix = args[iarg+1]
iarg += 2
else: error()
@ -51,32 +58,79 @@ while iarg < nargs:
cwd = os.getcwd()
lib = os.path.basename(cwd)
# create Makefile.auto as copy of Makefile.machine
# reset EXTRAMAKE if requested
def get_lammps_machine_flags(machine):
"""Parse Makefile.machine from LAMMPS, return dictionary of compiler flags"""
if not os.path.exists("../../src/MAKE/MACHINES/Makefile.%s" % machine):
error("Cannot locate src/MAKE/MACHINES/Makefile.%s" % machine)
lines = open("../../src/MAKE/MACHINES/Makefile.%s" % machine,
'r').readlines()
machine_flags = {}
for line in lines:
line = line.partition('#')[0]
line = line.rstrip()
words = line.split()
if (len(words) > 2):
if ((words[0] == 'CC') or (words[0] == 'CCFLAGS') or
(words[0] == 'SHFLAGS') or (words[0] == 'ARCHIVE') or
(words[0] == 'ARFLAGS') or (words[0] == 'SHELL')):
machine_flags[words[0]] = ' '.join(words[2:])
return machine_flags
def gen_colvars_makefile_machine(machine, machine_flags):
"""Generate Makefile.machine for Colvars given the compiler flags"""
machine_makefile = open("Makefile.%s" % machine, 'w')
machine_makefile.write('''# -*- makefile -*- to build Colvars module with %s
COLVARS_LIB = libcolvars.a
COLVARS_OBJ_DIR =
CXX = %s
CXXFLAGS = %s %s
AR = %s
ARFLAGS = %s
SHELL = %s
include Makefile.common
.PHONY: default clean
default: $(COLVARS_LIB) Makefile.lammps
clean:
-rm -f $(COLVARS_OBJS) $(COLVARS_LIB)
''' % (machine, machine_flags['CC'],
machine_flags['CCFLAGS'], machine_flags['SHFLAGS'] ,
machine_flags['ARCHIVE'], machine_flags['ARFLAGS'],
machine_flags['SHELL']))
if not os.path.exists("Makefile.%s" % machine):
machine_flags = get_lammps_machine_flags(machine)
gen_colvars_makefile_machine(machine, machine_flags)
if not os.path.exists("Makefile.%s" % machine):
error("lib/%s/Makefile.%s does not exist" % (lib,machine))
# create Makefile.auto as copy of Makefile.machine
# reset EXTRAMAKE if requested
lines = open("Makefile.%s" % machine,'r').readlines()
fp = open("Makefile.auto",'w')
for line in lines:
words = line.split()
if len(words) == 3 and extraflag and \
words[0] == "EXTRAMAKE" and words[1] == '=':
line = line.replace(words[2],"Makefile.lammps.%s" % suffix)
print >>fp,line,
fp.write(line)
fp.close()
# make the library via Makefile.auto
print "Building lib%s.a ..." % lib
cmd = "make -f Makefile.auto clean; make -f Makefile.auto"
txt = commands.getoutput(cmd)
print txt
print("Building lib%s.a ..." % lib)
cmd = ["make -f Makefile.auto clean"]
print(subprocess.check_output(cmd, shell=True))
cmd = ["make -f Makefile.auto -j12"]
print(subprocess.check_output(cmd, shell=True))
if os.path.exists("lib%s.a" % lib): print "Build was successful"
if os.path.exists("lib%s.a" % lib): print("Build was successful")
else: error("Build of lib/%s/lib%s.a was NOT successful" % (lib,lib))
if not os.path.exists("Makefile.lammps"):
print "lib/%s/Makefile.lammps was NOT created" % lib
print("lib/%s/Makefile.lammps was NOT created" % lib)

119
lib/colvars/Makefile.colvars
View File

@ -1,119 +0,0 @@
# library build -*- makefile -*- for colvars module
# which file will be copied to Makefile.lammps
EXTRAMAKE = Makefile.lammps.empty
# ------ SETTINGS ------
CXX = g++
CXXFLAGS = -O2 -g -Wall -fPIC -funroll-loops # -DCOLVARS_DEBUG
ARCHIVE = ar
ARCHFLAG = -rscv
SHELL = /bin/sh
# ------ DEFINITIONS ------
SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \
colvarbias_histogram.cpp colvarbias_meta.cpp colvarbias_restraint.cpp \
colvarcomp_angles.cpp colvarcomp_coordnums.cpp colvarcomp.cpp \
colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \
colvardeps.cpp colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \
colvarscript.cpp colvartypes.cpp colvarvalue.cpp
LIB = libcolvars.a
OBJ = $(SRC:.cpp=.o)
EXE = #colvars_standalone
# ------ MAKE PROCEDURE ------
default: $(LIB) $(EXE) Makefile.lammps
Makefile.lammps:
@cp $(EXTRAMAKE) Makefile.lammps
$(LIB): $(OBJ)
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB)
$(CXX) -o $@ $(CXXFLAGS) $^
# ------ MAKE FLAGS ------
.SUFFIXES:
.SUFFIXES: .cpp .o
.PHONY: default clean
# ------ COMPILE RULES ------
.cpp.o:
$(CXX) $(CXXFLAGS) -c $<
# ------ DEPENDENCIES ------
#
colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvaratoms.h
colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvar.h colvarparse.h colvardeps.h \
colvarbias_abf.h colvarbias.h colvargrid.h
colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarbias_alb.h colvar.h colvarparse.h \
colvardeps.h colvarbias_restraint.h colvarbias.h
colvarbias.o: colvarbias.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarbias.h colvar.h colvarparse.h colvardeps.h
colvarbias_histogram.o: colvarbias_histogram.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvar.h colvarparse.h \
colvardeps.h colvarbias_histogram.h colvarbias.h colvargrid.h
colvarbias_meta.o: colvarbias_meta.cpp colvar.h colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvarbias_meta.h colvarbias.h colvargrid.h
colvarbias_restraint.o: colvarbias_restraint.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarbias_restraint.h \
colvarbias.h colvar.h colvarparse.h colvardeps.h
colvarcomp_angles.o: colvarcomp_angles.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvar.h colvarparse.h colvardeps.h \
colvarcomp.h colvaratoms.h
colvarcomp_coordnums.o: colvarcomp_coordnums.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvaratoms.h colvar.h colvarcomp.h
colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvar.h colvarparse.h colvardeps.h colvarcomp.h \
colvaratoms.h
colvarcomp_distances.o: colvarcomp_distances.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvar.h colvarcomp.h colvaratoms.h
colvarcomp_protein.o: colvarcomp_protein.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h colvar.h \
colvarcomp.h colvaratoms.h
colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvar.h colvarcomp.h colvaratoms.h
colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
colvaratoms.h colvarscript.h colvarbias.h
colvardeps.o: colvardeps.cpp colvardeps.h colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarparse.h
colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
colvaratoms.h colvargrid.h
colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h colvar.h \
colvarbias.h colvarbias_abf.h colvargrid.h colvarbias_alb.h \
colvarbias_restraint.h colvarbias_histogram.h colvarbias_meta.h \
colvarscript.h
colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h
colvarscript.o: colvarscript.cpp colvarscript.h colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarbias.h colvar.h \
colvarparse.h colvardeps.h
colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h
colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h
# ------ CLEAN ------
clean:
-rm *.o *~ $(LIB)

65
lib/colvars/Makefile.common Normal file
View File

@ -0,0 +1,65 @@
# Shared -*- makefile -*- for multiple architectures
# Detect settings from PYTHON package (if defined)
include ../../src/Makefile.package.settings
ifeq ($(python_SYSINC),)
COLVARS_PYTHON_INCFLAGS =
else
COLVARS_PYTHON_INCFLAGS = -DCOLVARS_PYTHON $(python_SYSINC)
endif
# Detect debug settings
ifeq ($(COLVARS_DEBUG),)
COLVARS_DEBUG_INCFLAGS =
else
COLVARS_DEBUG_INCFLAGS= -DCOLVARS_DEBUG
endif
COLVARS_INCFLAGS = $(COLVARS_DEBUG_INCFLAGS) $(COLVARS_PYTHON_INCFLAGS)
.SUFFIXES:
.SUFFIXES: .cpp .o
COLVARS_SRCS = \
colvaratoms.cpp \
colvarbias_abf.cpp \
colvarbias_alb.cpp \
colvarbias.cpp \
colvarbias_histogram.cpp \
colvarbias_meta.cpp \
colvarbias_restraint.cpp \
colvarcomp_angles.cpp \
colvarcomp_coordnums.cpp \
colvarcomp.cpp \
colvarcomp_distances.cpp \
colvarcomp_protein.cpp \
colvarcomp_rotations.cpp \
colvar.cpp \
colvardeps.cpp \
colvargrid.cpp \
colvarmodule.cpp \
colvarparse.cpp \
colvarproxy.cpp \
colvarscript.cpp \
colvartypes.cpp \
colvarvalue.cpp
COLVARS_OBJS = $(COLVARS_SRCS:.cpp=.o)
.cpp.o:
$(CXX) $(CXXFLAGS) $(COLVARS_INCFLAGS) -c $<
$(COLVARS_LIB): Makefile.deps $(COLVARS_OBJS)
$(AR) $(ARFLAGS) $(COLVARS_LIB) $(COLVARS_OBJS)
Makefile.deps: $(COLVARS_SRCS)
@echo > $@
@for src in $^ ; do \
obj=`basename $$src .cpp`.o ; \
$(CXX) -MM $(COLVARS_INCFLAGS) \
-MT '$$(COLVARS_OBJ_DIR)'$$obj $$src >> $@ ; \
done
include Makefile.deps

78
lib/colvars/Makefile.deps Normal file
View File

@ -0,0 +1,78 @@
$(COLVARS_OBJ_DIR)colvaratoms.o: colvaratoms.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h \
colvarparse.h colvaratoms.h colvardeps.h
$(COLVARS_OBJ_DIR)colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h colvar.h \
colvarparse.h colvardeps.h colvarbias_abf.h colvarbias.h colvargrid.h
$(COLVARS_OBJ_DIR)colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h \
colvarbias_alb.h colvar.h colvarparse.h colvardeps.h colvarbias.h
$(COLVARS_OBJ_DIR)colvarbias.o: colvarbias.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h colvarbias.h \
colvar.h colvarparse.h colvardeps.h
$(COLVARS_OBJ_DIR)colvarbias_histogram.o: colvarbias_histogram.cpp \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvar.h colvarparse.h colvardeps.h colvarbias_histogram.h \
colvarbias.h colvargrid.h
$(COLVARS_OBJ_DIR)colvarbias_meta.o: colvarbias_meta.cpp colvar.h \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvarbias_meta.h colvarbias.h \
colvargrid.h
$(COLVARS_OBJ_DIR)colvarbias_restraint.o: colvarbias_restraint.cpp \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarbias_restraint.h colvarbias.h colvar.h colvarparse.h \
colvardeps.h
$(COLVARS_OBJ_DIR)colvarcomp_angles.o: colvarcomp_angles.cpp \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvar.h colvarparse.h colvardeps.h colvarcomp.h \
colvaratoms.h
$(COLVARS_OBJ_DIR)colvarcomp_coordnums.o: colvarcomp_coordnums.cpp \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvaratoms.h colvardeps.h colvar.h \
colvarcomp.h
$(COLVARS_OBJ_DIR)colvarcomp.o: colvarcomp.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h colvar.h \
colvarparse.h colvardeps.h colvarcomp.h colvaratoms.h
$(COLVARS_OBJ_DIR)colvarcomp_distances.o: colvarcomp_distances.cpp \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvar.h colvardeps.h colvarcomp.h \
colvaratoms.h
$(COLVARS_OBJ_DIR)colvarcomp_protein.o: colvarcomp_protein.cpp \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvar.h colvardeps.h colvarcomp.h \
colvaratoms.h
$(COLVARS_OBJ_DIR)colvarcomp_rotations.o: colvarcomp_rotations.cpp \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvar.h colvardeps.h colvarcomp.h \
colvaratoms.h
$(COLVARS_OBJ_DIR)colvar.o: colvar.cpp colvarmodule.h colvars_version.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvar.h \
colvardeps.h colvarcomp.h colvaratoms.h colvarscript.h colvarbias.h
$(COLVARS_OBJ_DIR)colvardeps.o: colvardeps.cpp colvardeps.h \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h
$(COLVARS_OBJ_DIR)colvargrid.o: colvargrid.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h \
colvarparse.h colvar.h colvardeps.h colvarcomp.h colvaratoms.h \
colvargrid.h
$(COLVARS_OBJ_DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h \
colvarparse.h colvar.h colvardeps.h colvarbias.h colvarbias_abf.h \
colvargrid.h colvarbias_alb.h colvarbias_histogram.h colvarbias_meta.h \
colvarbias_restraint.h colvarscript.h colvaratoms.h
$(COLVARS_OBJ_DIR)colvarparse.o: colvarparse.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h \
colvarparse.h
$(COLVARS_OBJ_DIR)colvarproxy.o: colvarproxy.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h \
colvarscript.h colvarbias.h colvar.h colvarparse.h colvardeps.h \
colvaratoms.h
$(COLVARS_OBJ_DIR)colvarscript.o: colvarscript.cpp colvarscript.h \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarbias.h colvar.h colvarparse.h colvardeps.h
$(COLVARS_OBJ_DIR)colvartypes.o: colvartypes.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h \
colvarparse.h
$(COLVARS_OBJ_DIR)colvarvalue.o: colvarvalue.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h

120
lib/colvars/Makefile.fermi
View File

@ -1,120 +0,0 @@
# library build -*- makefile -*- for colvars module
# which file will be copied to Makefile.lammps
EXTRAMAKE = Makefile.lammps.empty
# ------ SETTINGS ------
CXX = g++
CXXFLAGS = -O2 -mpc64 -g -fPIC \
-Wall -Wno-sign-compare # -DCOLVARS_DEBUG
ARCHIVE = ar
ARCHFLAG = -rscv
SHELL = /bin/sh
# ------ DEFINITIONS ------
SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \
colvarbias_histogram.cpp colvarbias_meta.cpp colvarbias_restraint.cpp \
colvarcomp_angles.cpp colvarcomp_coordnums.cpp colvarcomp.cpp \
colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \
colvardeps.cpp colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \
colvarscript.cpp colvartypes.cpp colvarvalue.cpp
LIB = libcolvars.a
OBJ = $(SRC:.cpp=.o)
EXE = #colvars_standalone
# ------ MAKE PROCEDURE ------
default: $(LIB) $(EXE) Makefile.lammps
Makefile.lammps:
@cp $(EXTRAMAKE) Makefile.lammps
$(LIB): $(OBJ)
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB)
$(CXX) -o $@ $(CXXFLAGS) $^
# ------ MAKE FLAGS ------
.SUFFIXES:
.SUFFIXES: .cpp .o
.PHONY: default clean
# ------ COMPILE RULES ------
.cpp.o:
$(CXX) $(CXXFLAGS) -c $<
# ------ DEPENDENCIES ------
#
colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvaratoms.h
colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvar.h colvarparse.h colvardeps.h \
colvarbias_abf.h colvarbias.h colvargrid.h
colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarbias_alb.h colvar.h colvarparse.h \
colvardeps.h colvarbias_restraint.h colvarbias.h
colvarbias.o: colvarbias.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarbias.h colvar.h colvarparse.h colvardeps.h
colvarbias_histogram.o: colvarbias_histogram.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvar.h colvarparse.h \
colvardeps.h colvarbias_histogram.h colvarbias.h colvargrid.h
colvarbias_meta.o: colvarbias_meta.cpp colvar.h colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvarbias_meta.h colvarbias.h colvargrid.h
colvarbias_restraint.o: colvarbias_restraint.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarbias_restraint.h \
colvarbias.h colvar.h colvarparse.h colvardeps.h
colvarcomp_angles.o: colvarcomp_angles.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvar.h colvarparse.h colvardeps.h \
colvarcomp.h colvaratoms.h
colvarcomp_coordnums.o: colvarcomp_coordnums.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvaratoms.h colvar.h colvarcomp.h
colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvar.h colvarparse.h colvardeps.h colvarcomp.h \
colvaratoms.h
colvarcomp_distances.o: colvarcomp_distances.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvar.h colvarcomp.h colvaratoms.h
colvarcomp_protein.o: colvarcomp_protein.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h colvar.h \
colvarcomp.h colvaratoms.h
colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvar.h colvarcomp.h colvaratoms.h
colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
colvaratoms.h colvarscript.h colvarbias.h
colvardeps.o: colvardeps.cpp colvardeps.h colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarparse.h
colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
colvaratoms.h colvargrid.h
colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h colvar.h \
colvarbias.h colvarbias_abf.h colvargrid.h colvarbias_alb.h \
colvarbias_restraint.h colvarbias_histogram.h colvarbias_meta.h \
colvarscript.h
colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h
colvarscript.o: colvarscript.cpp colvarscript.h colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarbias.h colvar.h \
colvarparse.h colvardeps.h
colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h
colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h
# ------ CLEAN ------
clean:
-rm *.o *~ $(LIB)

118
lib/colvars/Makefile.g++
View File

@ -1,119 +1,25 @@
# library build -*- makefile -*- for colvars module
# which file will be copied to Makefile.lammps
# -*- makefile -*- to build Colvars module with GNU compiler
EXTRAMAKE = Makefile.lammps.empty
# ------ SETTINGS ------
COLVARS_LIB = libcolvars.a
COLVARS_OBJ_DIR =
CXX = g++
CXXFLAGS = -O2 -g -fPIC -funroll-loops # -DCOLVARS_DEBUG
ARCHIVE = ar
ARCHFLAG = -rscv
CXXFLAGS = -O2 -g -Wall -fPIC -funroll-loops
AR = ar
ARFLAGS = -rscv
SHELL = /bin/sh
# ------ DEFINITIONS ------
SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \
colvarbias_histogram.cpp colvarbias_meta.cpp colvarbias_restraint.cpp \
colvarcomp_angles.cpp colvarcomp_coordnums.cpp colvarcomp.cpp \
colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \
colvardeps.cpp colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \
colvarscript.cpp colvartypes.cpp colvarvalue.cpp
LIB = libcolvars.a
OBJ = $(SRC:.cpp=.o)
EXE = #colvars_standalone
# ------ MAKE PROCEDURE ------
default: $(LIB) $(EXE) Makefile.lammps
Makefile.lammps:
@cp $(EXTRAMAKE) Makefile.lammps
$(LIB): $(OBJ)
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB)
$(CXX) -o $@ $(CXXFLAGS) $^
# ------ MAKE FLAGS ------
.SUFFIXES:
.SUFFIXES: .cpp .o
include Makefile.common
.PHONY: default clean
# ------ COMPILE RULES ------
.cpp.o:
$(CXX) $(CXXFLAGS) -c $<
# ------ DEPENDENCIES ------
#
colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvaratoms.h
colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvar.h colvarparse.h colvardeps.h \
colvarbias_abf.h colvarbias.h colvargrid.h
colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarbias_alb.h colvar.h colvarparse.h \
colvardeps.h colvarbias_restraint.h colvarbias.h
colvarbias.o: colvarbias.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarbias.h colvar.h colvarparse.h colvardeps.h
colvarbias_histogram.o: colvarbias_histogram.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvar.h colvarparse.h \
colvardeps.h colvarbias_histogram.h colvarbias.h colvargrid.h
colvarbias_meta.o: colvarbias_meta.cpp colvar.h colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvarbias_meta.h colvarbias.h colvargrid.h
colvarbias_restraint.o: colvarbias_restraint.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarbias_restraint.h \
colvarbias.h colvar.h colvarparse.h colvardeps.h
colvarcomp_angles.o: colvarcomp_angles.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvar.h colvarparse.h colvardeps.h \
colvarcomp.h colvaratoms.h
colvarcomp_coordnums.o: colvarcomp_coordnums.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvaratoms.h colvar.h colvarcomp.h
colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvar.h colvarparse.h colvardeps.h colvarcomp.h \
colvaratoms.h
colvarcomp_distances.o: colvarcomp_distances.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvar.h colvarcomp.h colvaratoms.h
colvarcomp_protein.o: colvarcomp_protein.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h colvar.h \
colvarcomp.h colvaratoms.h
colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvar.h colvarcomp.h colvaratoms.h
colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
colvaratoms.h colvarscript.h colvarbias.h
colvardeps.o: colvardeps.cpp colvardeps.h colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarparse.h
colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
colvaratoms.h colvargrid.h
colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h colvar.h \
colvarbias.h colvarbias_abf.h colvargrid.h colvarbias_alb.h \
colvarbias_restraint.h colvarbias_histogram.h colvarbias_meta.h \
colvarscript.h
colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h
colvarscript.o: colvarscript.cpp colvarscript.h colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarbias.h colvar.h \
colvarparse.h colvardeps.h
colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h
colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h
# ------ CLEAN ------
default: $(COLVARS_LIB) Makefile.lammps
clean:
-rm *.o *~ $(LIB)
-rm -f $(COLVARS_OBJS) $(COLVARS_LIB)
Makefile.lammps:
-cp $(EXTRAMAKE) Makefile.lammps

5
lib/colvars/Makefile.g++-debug Normal file
View File

@ -0,0 +1,5 @@
# -*- makefile -*- to build Colvars module with GNU compiler
COLVARS_DEBUG = "YES"
include Makefile.g++

5
lib/colvars/Makefile.lammps Normal file
View File

@ -0,0 +1,5 @@
# Settings that the LAMMPS build will import when this package library is used
colvars_SYSINC =
colvars_SYSLIB =
colvars_SYSPATH =

2
lib/colvars/Makefile.lammps.debug
View File

@ -1,5 +1,5 @@
# Settings that the LAMMPS build will import when this package library is used
colvars_SYSINC = # -DCOLVARS_DEBUG
colvars_SYSINC = -DCOLVARS_DEBUG
colvars_SYSLIB =
colvars_SYSPATH =

2
lib/colvars/Makefile.lammps.empty
View File

@ -1,5 +1,5 @@
# Settings that the LAMMPS build will import when this package library is used
colvars_SYSINC = # -DCOLVARS_DEBUG
colvars_SYSINC =
colvars_SYSLIB =
colvars_SYSPATH =

124
lib/colvars/Makefile.mingw32-cross
View File

@ -1,127 +1,31 @@
# library build -*- makefile -*- for colvars module
# which file will be copied to Makefile.lammps
# -*- makefile -*- to build Colvars module with MinGW 32-bit
EXTRAMAKE = Makefile.lammps.empty
# ------ SETTINGS ------
COLVARS_LIB = libcolvars.a
COLVARS_OBJ_DIR = Obj_mingw64/
CXX = i686-w64-mingw32-g++
CXXFLAGS = -O2 -march=i686 -mtune=generic -mfpmath=387 -mpc64 \
-fno-rtti -fno-exceptions -finline-functions \
-ffast-math -funroll-loops -fstrict-aliasing \
-Wall -W -Wno-uninitialized
ARCHIVE = i686-w64-mingw32-ar
ARCHFLAG = -rscv
AR = i686-w64-mingw32-ar
ARFLAGS = -rscv
SHELL = /bin/sh
# ------ DEFINITIONS ------
SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \
colvarbias_histogram.cpp colvarbias_meta.cpp colvarbias_restraint.cpp \
colvarcomp_angles.cpp colvarcomp_coordnums.cpp colvarcomp.cpp \
colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \
colvardeps.cpp colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \
colvarscript.cpp colvartypes.cpp colvarvalue.cpp
DIR = Obj_mingw32/
LIB = $(DIR)libcolvars.a
OBJ = $(SRC:%.cpp=$(DIR)%.o)
EXE = #colvars_standalone
# ------ MAKE PROCEDURE ------
default: $(DIR) $(LIB) $(EXE) Makefile.lammps
$(DIR):
mkdir $(DIR)
Makefile.lammps:
@cp $(EXTRAMAKE) Makefile.lammps
$(LIB): $(DIR) $(OBJ)
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
@cp $(EXTRAMAKE) Makefile.lammps
$(DIR)colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB)
$(CXX) -o $@ $(CXXFLAGS) $^
# ------ MAKE FLAGS ------
.SUFFIXES:
.SUFFIXES: .cpp .o
include Makefile.common
.PHONY: default clean
# ------ COMPILE RULES ------
default: $(COLVARS_OBJ_DIR) $(COLVARS_LIB) Makefile.lammps
$(DIR)%.o: %.cpp
$(CXX) $(CXXFLAGS) -c $< -o $@
# ------ DEPENDENCIES ------
#
$(DIR)colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvaratoms.h
$(DIR)colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvar.h colvarparse.h colvardeps.h \
colvarbias_abf.h colvarbias.h colvargrid.h
$(DIR)colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarbias_alb.h colvar.h colvarparse.h \
colvardeps.h colvarbias_restraint.h colvarbias.h
$(DIR)colvarbias.o: colvarbias.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarbias.h colvar.h colvarparse.h colvardeps.h
$(DIR)colvarbias_histogram.o: colvarbias_histogram.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvar.h colvarparse.h \
colvardeps.h colvarbias_histogram.h colvarbias.h colvargrid.h
$(DIR)colvarbias_meta.o: colvarbias_meta.cpp colvar.h colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvarbias_meta.h colvarbias.h colvargrid.h
$(DIR)colvarbias_restraint.o: colvarbias_restraint.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarbias_restraint.h \
colvarbias.h colvar.h colvarparse.h colvardeps.h
$(DIR)colvarcomp_angles.o: colvarcomp_angles.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvar.h colvarparse.h colvardeps.h \
colvarcomp.h colvaratoms.h
$(DIR)colvarcomp_coordnums.o: colvarcomp_coordnums.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvaratoms.h colvar.h colvarcomp.h
$(DIR)colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvar.h colvarparse.h colvardeps.h colvarcomp.h \
colvaratoms.h
$(DIR)colvarcomp_distances.o: colvarcomp_distances.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvar.h colvarcomp.h colvaratoms.h
$(DIR)colvarcomp_protein.o: colvarcomp_protein.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h colvar.h \
colvarcomp.h colvaratoms.h
$(DIR)colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvar.h colvarcomp.h colvaratoms.h
$(DIR)colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
colvaratoms.h colvarscript.h colvarbias.h
$(DIR)colvardeps.o: colvardeps.cpp colvardeps.h colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarparse.h
$(DIR)colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
colvaratoms.h colvargrid.h
$(DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h colvar.h \
colvarbias.h colvarbias_abf.h colvargrid.h colvarbias_alb.h \
colvarbias_restraint.h colvarbias_histogram.h colvarbias_meta.h \
colvarscript.h
$(DIR)colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h
$(DIR)colvarscript.o: colvarscript.cpp colvarscript.h colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarbias.h colvar.h \
colvarparse.h colvardeps.h
$(DIR)colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h
$(DIR)colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h
# ------ CLEAN ------
$(COLVARS_OBJ_DIR):
mkdir $(COLVARS_OBJ_DIR)
clean:
-rm $(DIR)*.o *~ $(LIB)
-rmdir $(DIR)
-rm -f $(COLVARS_OBJS) $(COLVARS_LIB)
-rmdir $(COLVARS_OBJ_DIR)
Makefile.lammps:
-cp $(EXTRAMAKE) Makefile.lammps

124
lib/colvars/Makefile.mingw64-cross
View File

@ -1,127 +1,31 @@
# library build -*- makefile -*- for colvars module
# which file will be copied to Makefile.lammps
# -*- makefile -*- to build Colvars module with MinGW 32-bit
EXTRAMAKE = Makefile.lammps.empty
# ------ SETTINGS ------
COLVARS_LIB = libcolvars.a
COLVARS_OBJ_DIR = Obj_mingw32/
CXX = x86_64-w64-mingw32-g++
CXXFLAGS = -O2 -march=core2 -mtune=core2 -mpc64 -msse2 \
-fno-rtti -fno-exceptions -finline-functions \
-ffast-math -funroll-loops -fstrict-aliasing \
-Wall -W -Wno-uninitialized
ARCHIVE = x86_64-w64-mingw32-ar
ARCHFLAG = -rscv
AR = x86_64-w64-mingw32-ar
ARFLAGS = -rscv
SHELL = /bin/sh
# ------ DEFINITIONS ------
SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \
colvarbias_histogram.cpp colvarbias_meta.cpp colvarbias_restraint.cpp \
colvarcomp_angles.cpp colvarcomp_coordnums.cpp colvarcomp.cpp \
colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \
colvardeps.cpp colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \
colvarscript.cpp colvartypes.cpp colvarvalue.cpp
DIR = Obj_mingw64/
LIB = $(DIR)libcolvars.a
OBJ = $(SRC:%.cpp=$(DIR)%.o)
EXE = #colvars_standalone
# ------ MAKE PROCEDURE ------
default: $(DIR) $(LIB) $(EXE) Makefile.lammps
$(DIR):
mkdir $(DIR)
Makefile.lammps:
@cp $(EXTRAMAKE) Makefile.lammps
$(LIB): $(DIR) $(OBJ)
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
@cp $(EXTRAMAKE) Makefile.lammps
$(DIR)colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB)
$(CXX) -o $@ $(CXXFLAGS) $^
# ------ MAKE FLAGS ------
.SUFFIXES:
.SUFFIXES: .cpp .o
include Makefile.common
.PHONY: default clean
# ------ COMPILE RULES ------
default: $(COLVARS_OBJ_DIR) $(COLVARS_LIB) Makefile.lammps
$(DIR)%.o: %.cpp
$(CXX) $(CXXFLAGS) -c $< -o $@
# ------ DEPENDENCIES ------
#
$(DIR)colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvaratoms.h
$(DIR)colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvar.h colvarparse.h colvardeps.h \
colvarbias_abf.h colvarbias.h colvargrid.h
$(DIR)colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarbias_alb.h colvar.h colvarparse.h \
colvardeps.h colvarbias_restraint.h colvarbias.h
$(DIR)colvarbias.o: colvarbias.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarbias.h colvar.h colvarparse.h colvardeps.h
$(DIR)colvarbias_histogram.o: colvarbias_histogram.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvar.h colvarparse.h \
colvardeps.h colvarbias_histogram.h colvarbias.h colvargrid.h
$(DIR)colvarbias_meta.o: colvarbias_meta.cpp colvar.h colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvarbias_meta.h colvarbias.h colvargrid.h
$(DIR)colvarbias_restraint.o: colvarbias_restraint.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarbias_restraint.h \
colvarbias.h colvar.h colvarparse.h colvardeps.h
$(DIR)colvarcomp_angles.o: colvarcomp_angles.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvar.h colvarparse.h colvardeps.h \
colvarcomp.h colvaratoms.h
$(DIR)colvarcomp_coordnums.o: colvarcomp_coordnums.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvaratoms.h colvar.h colvarcomp.h
$(DIR)colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvar.h colvarparse.h colvardeps.h colvarcomp.h \
colvaratoms.h
$(DIR)colvarcomp_distances.o: colvarcomp_distances.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvar.h colvarcomp.h colvaratoms.h
$(DIR)colvarcomp_protein.o: colvarcomp_protein.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h colvar.h \
colvarcomp.h colvaratoms.h
$(DIR)colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h \
colvar.h colvarcomp.h colvaratoms.h
$(DIR)colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
colvaratoms.h colvarscript.h colvarbias.h
$(DIR)colvardeps.o: colvardeps.cpp colvardeps.h colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarparse.h
$(DIR)colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
colvaratoms.h colvargrid.h
$(DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h \
colvarproxy.h colvarvalue.h colvarparse.h colvardeps.h colvar.h \
colvarbias.h colvarbias_abf.h colvargrid.h colvarbias_alb.h \
colvarbias_restraint.h colvarbias_histogram.h colvarbias_meta.h \
colvarscript.h
$(DIR)colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h
$(DIR)colvarscript.o: colvarscript.cpp colvarscript.h colvarmodule.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarbias.h colvar.h \
colvarparse.h colvardeps.h
$(DIR)colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h
$(DIR)colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \
colvarvalue.h
# ------ CLEAN ------
$(COLVARS_OBJ_DIR):
mkdir $(COLVARS_OBJ_DIR)
clean:
-rm $(DIR)*.o *~ $(LIB)
-rmdir $(DIR)
-rm -f $(COLVARS_OBJS) $(COLVARS_LIB)
-rmdir $(COLVARS_OBJ_DIR)
Makefile.lammps:
-cp $(EXTRAMAKE) Makefile.lammps

105
lib/colvars/README
View File

@ -1,49 +1,35 @@
This library is the portable "colvars" module, originally interfaced
with the NAMD MD code, to provide an extensible software framework,
that allows enhanced sampling in molecular dynamics simulations.
The module is written to maximize performance, portability,
flexibility of usage for the user, and extensibility for the developer.
## Collective variables module (Colvars)
The development of the colvars library is now hosted on github at:
http://colvars.github.io/
You can use this site to get access to the latest development sources
and the up-to-date documentation.
A software module for molecular simulation and analysis that provides a
high-performance implementation of sampling algorithms defined on a reduced
space of continuously differentiable functions (aka collective variables).
Copy of the specific documentation is also in
doc/PDF/colvars-refman-lammps.pdf
The module itself implements a variety of functions and algorithms, including
free-energy estimators based on thermodynamic forces, non-equilibrium work and
probability distributions.
Please report bugs and request new features at:
https://github.com/colvars/colvars/issues
For a brief description see:
http://colvars.github.io/
https://github.com/colvars/colvars/
The following publications describe the principles of
the implementation of this library:
Using collective variables to drive molecular dynamics simulations,
Giacomo Fiorin , Michael L. Klein & Jérôme Hénin (2013):
Molecular Physics DOI:10.1080/00268976.2013.813594
Exploring Multidimensional Free Energy Landscapes Using
Time-Dependent Biases on Collective Variables,
J. Hénin, G. Fiorin, C. Chipot, and M. L. Klein,
J. Chem. Theory Comput., 6, 35-47 (2010).
-------------------------------------------------
## How to build
This directory has source files to build a library that LAMMPS
links against when using the USER-COLVARS package.
This library must be built with a C++ compiler, before LAMMPS is
built, so LAMMPS can link against it.
This library must be built with a C++ compiler, *before* LAMMPS is built, so
that LAMMPS can link against it. You can use the provided Makefile.* files or
create your own, specific to your compiler and system. For example:
You can type "make lib-colvars" from the src directory to see help on
how to build this library via make commands, or you can do the same
thing by typing "python Install.py" from within this directory, or you
can do it manually by following the instructions below.
make -f Makefile.g++
Build the library using one of the provided Makefile.* files or create
your own, specific to your compiler and system. For example:
will use the GNU C++ compiler and is a good template to start.
make -f Makefile.g++
**Optional**: if you use the Install.py script provided in this folder, you
can give the machine name as the '-m' argument. This can be the suffix of one
of the files from either this folder, or from src/MAKE.
*This is only supported by the Install.py within the lib/colvars folder*.
When you are done building this library, two files should
exist in this directory:
@ -51,23 +37,42 @@ exist in this directory:
libcolvars.a the library LAMMPS will link against
Makefile.lammps settings the LAMMPS Makefile will import
Makefile.lammps is created by the make command, by copying one of the
Makefile.lammps.* files. See the EXTRAMAKE setting at the top of the
Makefile.* files.
IMPORTANT: You must examine the final Makefile.lammps to insure it is
correct for your system, else the LAMMPS build will likely fail.
Makefile.lammps has settings for 3 variables:
user-colvars_SYSINC = leave blank for this package unless debugging
user-colvars_SYSLIB = leave blank for this package
user-colvars_SYSPATH = leave blank for this package
You have several choices for these settings:
Since they do not normally need to be set, the settings in
Makefile.lammps.empty should work.
If you want to set a debug flag recognized by the library, the
settings in Makefile.lammps.debug should work.
settings in Makefile.common should work.
## Documentation
For the reference manual see:
http://colvars.github.io/colvars-refman-lammps
A copy of reference manual is also in:
doc/PDF/colvars-refman-lammps.pdf
Also included is a Doxygen-based developer documentation:
http://colvars.github.io/doxygen/html/
The reference article is:
G. Fiorin, M. L. Klein, and J. Henin,
Molecular Physics 111, 3345 (2013).
http://dx.doi.org/10.1080/00268976.2013.813594
## Updating to the latest version
To recompile LAMMPS with the most recent version of this module, the `master`
branch of this repository from GitHub, or clone it via git:
git clone https://github.com/colvars/colvars.git
and run the provided `update-colvars-code.sh` script against the unpacked
LAMMPS source tree:
./update-colvars-code.sh /path/to/lammps/folder
Please report bugs and request new features at:
https://github.com/colvars/colvars/issues

400
lib/colvars/colvar.cpp
View File

@ -1,5 +1,5 @@
// -*- c++ -*-
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/colvars/colvars
@ -7,13 +7,14 @@
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
#include "colvar.h"
#include "colvarcomp.h"
#include "colvarscript.h"
// used in build_atom_list()
#include <algorithm>
@ -25,8 +26,10 @@ bool compare(colvar::cvc *i, colvar::cvc *j) {
colvar::colvar()
: prev_timestep(-1)
{
// Initialize static array once and for all
runave_os = NULL;
init_cv_requires();
}
@ -66,6 +69,13 @@ int colvar::init(std::string const &conf)
size_t i;
#ifdef LEPTON
error_code |= init_custom_function(conf);
if (error_code != COLVARS_OK) {
return cvm::get_error();
}
#endif
// Setup colvar as scripted function of components
if (get_keyval(conf, "scriptedFunction", scripted_function,
"", colvarparse::parse_silent)) {
@ -122,7 +132,7 @@ int colvar::init(std::string const &conf)
}
}
if (!is_enabled(f_cv_scripted)) {
if (!(is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function))) {
colvarvalue const &cvc_value = (cvcs[0])->value();
if (cvm::debug())
cvm::log ("This collective variable is a "+
@ -141,7 +151,7 @@ int colvar::init(std::string const &conf)
// check for linear combinations
{
bool lin = !is_enabled(f_cv_scripted);
bool lin = !(is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function));
for (i = 0; i < cvcs.size(); i++) {
// FIXME this is a reverse dependency, ie. cv feature depends on cvc flag
@ -206,7 +216,7 @@ int colvar::init(std::string const &conf)
for (i = 0; i < cvcs.size(); i++) {
// components may have different types only for scripted functions
if (!is_enabled(f_cv_scripted) && (colvarvalue::check_types(cvcs[i]->value(),
if (!(is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function)) && (colvarvalue::check_types(cvcs[i]->value(),
cvcs[0]->value())) ) {
cvm::error("ERROR: you are definining this collective variable "
"by using components of different types. "
@ -223,7 +233,6 @@ int colvar::init(std::string const &conf)
// at this point, the colvar's type is defined
f.type(value());
f_accumulated.type(value());
x_old.type(value());
v_fdiff.type(value());
@ -239,18 +248,23 @@ int colvar::init(std::string const &conf)
reset_bias_force();
get_keyval(conf, "timeStepFactor", time_step_factor, 1);
if (time_step_factor < 0) {
cvm::error("Error: timeStepFactor must be positive.\n");
return COLVARS_ERROR;
}
if (time_step_factor != 1) {
enable(f_cv_multiple_ts);
}
// TODO use here information from the CVCs' own natural boundaries
error_code |= init_grid_parameters(conf);
get_keyval(conf, "timeStepFactor", time_step_factor, 1);
error_code |= init_extended_Lagrangian(conf);
error_code |= init_output_flags(conf);
// Start in active state by default
// Now that the children are defined we can solve dependencies
enable(f_cv_active);
// Make sure dependency side-effects are correct
refresh_deps();
if (cvm::b_analysis)
parse_analysis(conf);
@ -262,6 +276,158 @@ int colvar::init(std::string const &conf)
}
#ifdef LEPTON
int colvar::init_custom_function(std::string const &conf)
{
std::string expr;
std::vector<Lepton::ParsedExpression> pexprs;
Lepton::ParsedExpression pexpr;
size_t pos = 0; // current position in config string
double *ref;
if (!key_lookup(conf, "customFunction", &expr, &pos)) {
return COLVARS_OK;
}
enable(f_cv_custom_function);
cvm::log("This colvar uses a custom function.\n");
do {
if (cvm::debug())
cvm::log("Parsing expression \"" + expr + "\".\n");
try {
pexpr = Lepton::Parser::parse(expr);
pexprs.push_back(pexpr);
}
catch (...) {
cvm::error("Error parsing expression \"" + expr + "\".\n", INPUT_ERROR);
return INPUT_ERROR;
}
try {
value_evaluators.push_back(
new Lepton::CompiledExpression(pexpr.createCompiledExpression()));
// Define variables for cvc values
// Stored in order: expr1, cvc1, cvc2, expr2, cvc1...
for (size_t i = 0; i < cvcs.size(); i++) {
for (size_t j = 0; j < cvcs[i]->value().size(); j++) {
std::string vn = cvcs[i]->name +
(cvcs[i]->value().size() > 1 ? cvm::to_str(j+1) : "");
try {
ref =&value_evaluators.back()->getVariableReference(vn);
}
catch (...) { // Variable is absent from expression
// To keep the same workflow, we use a pointer to a double here
// that will receive CVC values - even though none was allocated by Lepton
ref = &dev_null;
if (cvm::debug())
cvm::log("Variable " + vn + " is absent from expression \"" + expr + "\".\n");
}
value_eval_var_refs.push_back(ref);
}
}
}
catch (...) {
cvm::error("Error compiling expression \"" + expr + "\".\n", INPUT_ERROR);
return INPUT_ERROR;
}
} while (key_lookup(conf, "customFunction", &expr, &pos));
// Now define derivative with respect to each scalar sub-component
for (size_t i = 0; i < cvcs.size(); i++) {
for (size_t j = 0; j < cvcs[i]->value().size(); j++) {
std::string vn = cvcs[i]->name +
(cvcs[i]->value().size() > 1 ? cvm::to_str(j+1) : "");
// Element ordering: we want the
// gradient vector of derivatives of all elements of the colvar
// wrt to a given element of a cvc ([i][j])
for (size_t c = 0; c < pexprs.size(); c++) {
gradient_evaluators.push_back(
new Lepton::CompiledExpression(pexprs[c].differentiate(vn).createCompiledExpression()));
// and record the refs to each variable in those expressions
for (size_t k = 0; k < cvcs.size(); k++) {
for (size_t l = 0; l < cvcs[k]->value().size(); l++) {
std::string vvn = cvcs[k]->name +
(cvcs[k]->value().size() > 1 ? cvm::to_str(l+1) : "");
try {
ref = &gradient_evaluators.back()->getVariableReference(vvn);
}
catch (...) { // Variable is absent from derivative
// To keep the same workflow, we use a pointer to a double here
// that will receive CVC values - even though none was allocated by Lepton
if (cvm::debug())
cvm::log("Variable " + vvn + " is absent from derivative of \"" + expr + "\" wrt " + vn + ".\n");
ref = &dev_null;
}
grad_eval_var_refs.push_back(ref);
}
}
}
}
}
if (value_evaluators.size() == 0) {
cvm::error("Error: no custom function defined.\n", INPUT_ERROR);
return INPUT_ERROR;
}
std::string type_str;
bool b_type_specified = get_keyval(conf, "customFunctionType",
type_str, "scalar", parse_silent);
x.type(colvarvalue::type_notset);
int t;
for (t = 0; t < colvarvalue::type_all; t++) {
if (type_str == colvarvalue::type_keyword(colvarvalue::Type(t))) {
x.type(colvarvalue::Type(t));
break;
}
}
if (x.type() == colvarvalue::type_notset) {
cvm::error("Could not parse custom colvar type.", INPUT_ERROR);
return INPUT_ERROR;
}
// Guess type based on number of expressions
if (!b_type_specified) {
if (value_evaluators.size() == 1) {
x.type(colvarvalue::type_scalar);
} else {
x.type(colvarvalue::type_vector);
}
}
if (x.type() == colvarvalue::type_vector) {
x.vector1d_value.resize(value_evaluators.size());
}
x_reported.type(x);
cvm::log(std::string("Expecting colvar value of type ")
+ colvarvalue::type_desc(x.type())
+ (x.type()==colvarvalue::type_vector ? " of size " + cvm::to_str(x.size()) : "")
+ ".\n");
if (x.size() != value_evaluators.size()) {
cvm::error("Error: based on custom function type, expected "
+ cvm::to_str(x.size()) + " scalar expressions, but "
+ cvm::to_str(value_evaluators.size() + " were found.\n"));
return INPUT_ERROR;
}
return COLVARS_OK;
}
#else
int colvar::init_custom_function(std::string const &conf)
{
return COLVARS_OK;
}
#endif // #ifdef LEPTON
int colvar::init_grid_parameters(std::string const &conf)
{
colvarmodule *cv = cvm::main();
@ -326,7 +492,8 @@ int colvar::init_grid_parameters(std::string const &conf)
std::string const walls_conf("\n\
harmonicWalls {\n\
name "+this->name+"w\n\
colvars "+this->name+"\n"+lw_conf+uw_conf+
colvars "+this->name+"\n"+lw_conf+uw_conf+"\
timeStepFactor "+cvm::to_str(time_step_factor)+"\n"+
"}\n");
cv->append_new_config(walls_conf);
}
@ -372,17 +539,14 @@ harmonicWalls {\n\
int colvar::init_extended_Lagrangian(std::string const &conf)
{
bool b_extended_Lagrangian;
get_keyval(conf, "extendedLagrangian", b_extended_Lagrangian, false);
get_keyval_feature(this, conf, "extendedLagrangian", f_cv_extended_Lagrangian, false);
if (b_extended_Lagrangian) {
if (is_enabled(f_cv_extended_Lagrangian)) {
cvm::real temp, tolerance, period;
cvm::log("Enabling the extended Lagrangian term for colvar \""+
this->name+"\".\n");
enable(f_cv_extended_Lagrangian);
xr.type(value());
vr.type(value());
fr.type(value());
@ -404,7 +568,7 @@ int colvar::init_extended_Lagrangian(std::string const &conf)
return INPUT_ERROR;
}
ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance);
cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " kcal/mol/U^2");
cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " [E]/U^2");
get_keyval(conf, "extendedTimeConstant", period, 200.0);
if (period <= 0.0) {
@ -412,7 +576,7 @@ int colvar::init_extended_Lagrangian(std::string const &conf)
}
ext_mass = (cvm::boltzmann() * temp * period * period)
/ (4.0 * PI * PI * tolerance * tolerance);
cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " kcal/mol/(U/fs)^2 (U: colvar unit)");
cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " [E]/(U/fs)^2 (U: colvar unit)");
{
bool b_output_energy;
@ -429,8 +593,9 @@ int colvar::init_extended_Lagrangian(std::string const &conf)
}
if (ext_gamma != 0.0) {
enable(f_cv_Langevin);
ext_gamma *= 1.0e-3; // convert from ps-1 to fs-1
ext_sigma = std::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / cvm::dt());
ext_gamma *= 1.0e-3; // correct as long as input is required in ps-1 and cvm::dt() is in fs
// Adjust Langevin sigma for slow time step if time_step_factor != 1
ext_sigma = std::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / (cvm::dt() * cvm::real(time_step_factor)));
}
}
@ -486,8 +651,8 @@ template<typename def_class_name> int colvar::init_components_type(std::string c
size_t pos = 0;
while ( this->key_lookup(conf,
def_config_key,
def_conf,
pos) ) {
&def_conf,
&pos) ) {
if (!def_conf.size()) continue;
cvm::log("Initializing "
"a new \""+std::string(def_config_key)+"\" component"+
@ -514,6 +679,7 @@ template<typename def_class_name> int colvar::init_components_type(std::string c
if ( (cvcp->period != 0.0) || (cvcp->wrap_center != 0.0) ) {
if ( (cvcp->function_type != std::string("distance_z")) &&
(cvcp->function_type != std::string("dihedral")) &&
(cvcp->function_type != std::string("polar_phi")) &&
(cvcp->function_type != std::string("spin_angle")) ) {
cvm::error("Error: invalid use of period and/or "
"wrapAround in a \""+
@ -566,6 +732,10 @@ int colvar::init_components(std::string const &conf)
"on an axis", "distanceZ");
error_code |= init_components_type<distance_xy>(conf, "distance projection "
"on a plane", "distanceXY");
error_code |= init_components_type<polar_theta>(conf, "spherical polar angle theta",
"polarTheta");
error_code |= init_components_type<polar_phi>(conf, "spherical azimuthal angle phi",
"polarPhi");
error_code |= init_components_type<distance_inv>(conf, "average distance "
"weighted by inverse power", "distanceInv");
error_code |= init_components_type<distance_pairs>(conf, "N1xN2-long vector "
@ -618,16 +788,18 @@ int colvar::init_components(std::string const &conf)
}
int colvar::refresh_deps()
void colvar::do_feature_side_effects(int id)
{
// If enabled features are changed upstream, the features below should be refreshed
if (is_enabled(f_cv_total_force_calc)) {
switch (id) {
case f_cv_total_force_calc:
cvm::request_total_force();
}
if (is_enabled(f_cv_collect_gradient) && atom_ids.size() == 0) {
break;
case f_cv_collect_gradient:
if (atom_ids.size() == 0) {
build_atom_list();
}
return COLVARS_OK;
break;
}
}
@ -688,15 +860,14 @@ int colvar::parse_analysis(std::string const &conf)
cvm::error("Error: runAveStride must be commensurate with the restart frequency.\n", INPUT_ERROR);
}
std::string runave_outfile;
get_keyval(conf, "runAveOutputFile", runave_outfile,
std::string(cvm::output_prefix()+"."+
this->name+".runave.traj"));
size_t const this_cv_width = x.output_width(cvm::cv_width);
cvm::backup_file(runave_outfile.c_str());
runave_os.open(runave_outfile.c_str());
runave_os << "# " << cvm::wrap_string("step", cvm::it_width-2)
cvm::proxy->backup_file(runave_outfile);
runave_os = cvm::proxy->output_stream(runave_outfile);
*runave_os << "# " << cvm::wrap_string("step", cvm::it_width-2)
<< " "
<< cvm::wrap_string("running average", this_cv_width)
<< " "
@ -768,6 +939,10 @@ void colvar::setup() {
colvar::~colvar()
{
// There is no need to call free_children_deps() here
// because the children are cvcs and will be deleted
// just below
// Clear references to this colvar's cvcs as children
// for dependency purposes
remove_all_children();
@ -792,6 +967,22 @@ colvar::~colvar()
break;
}
}
#ifdef LEPTON
for (std::vector<Lepton::CompiledExpression *>::iterator cei = value_evaluators.begin();
cei != value_evaluators.end();
++cei) {
if (*cei != NULL) delete (*cei);
}
value_evaluators.clear();
for (std::vector<Lepton::CompiledExpression *>::iterator gei = gradient_evaluators.begin();
gei != gradient_evaluators.end();
++gei) {
if (*gei != NULL) delete (*gei);
}
gradient_evaluators.clear();
#endif
}
@ -911,7 +1102,6 @@ int colvar::calc_cvc_values(int first_cvc, size_t num_cvcs)
int colvar::collect_cvc_values()
{
x.reset();
size_t i;
// combine them appropriately, using either a scripted function or a polynomial
if (is_enabled(f_cv_scripted)) {
@ -925,9 +1115,26 @@ int colvar::collect_cvc_values()
cvm::error("Error running scripted colvar");
return COLVARS_OK;
}
#ifdef LEPTON
} else if (is_enabled(f_cv_custom_function)) {
size_t l = 0; // index in the vector of variable references
for (size_t i = 0; i < x.size(); i++) {
// Fill Lepton evaluator variables with CVC values, serialized into scalars
for (size_t j = 0; j < cvcs.size(); j++) {
for (size_t k = 0; k < cvcs[j]->value().size(); k++) {
*(value_eval_var_refs[l++]) = cvcs[j]->value()[k];
}
}
x[i] = value_evaluators[i]->evaluate();
}
#endif
} else if (x.type() == colvarvalue::type_scalar) {
// polynomial combination allowed
for (i = 0; i < cvcs.size(); i++) {
for (size_t i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->is_enabled()) continue;
x += (cvcs[i])->sup_coeff *
( ((cvcs[i])->sup_np != 1) ?
@ -935,7 +1142,7 @@ int colvar::collect_cvc_values()
(cvcs[i])->value().real_value );
}
} else {
for (i = 0; i < cvcs.size(); i++) {
for (size_t i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->is_enabled()) continue;
x += (cvcs[i])->sup_coeff * (cvcs[i])->value();
}
@ -984,16 +1191,9 @@ int colvar::calc_cvc_gradients(int first_cvc, size_t num_cvcs)
(cvcs[i])->calc_gradients();
// if requested, propagate (via chain rule) the gradients above
// to the atoms used to define the roto-translation
// This could be integrated in the CVC base class
for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
if (cvcs[i]->atom_groups[ig]->b_fit_gradients)
cvcs[i]->atom_groups[ig]->calc_fit_gradients();
if (cvcs[i]->is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging gradients for " + cvcs[i]->description);
cvcs[i]->debug_gradients(cvcs[i]->atom_groups[ig]);
}
}
(cvcs[i])->calc_fit_gradients();
if ((cvcs[i])->is_enabled(f_cvc_debug_gradient))
(cvcs[i])->debug_gradients();
}
cvm::decrease_depth();
@ -1011,13 +1211,6 @@ int colvar::collect_cvc_gradients()
size_t i;
if (is_enabled(f_cv_collect_gradient)) {
if (is_enabled(f_cv_scripted)) {
cvm::error("Collecting atomic gradients is not implemented for "
"scripted colvars.", COLVARS_NOT_IMPLEMENTED);
return COLVARS_NOT_IMPLEMENTED;
}
// Collect the atomic gradients inside colvar object
for (unsigned int a = 0; a < atomic_gradients.size(); a++) {
atomic_gradients[a].reset();
@ -1214,6 +1407,11 @@ cvm::real colvar::update_forces_energy()
// set to zero the applied force
f.type(value());
f.reset();
fr.reset();
// If we are not active at this timestep, that's all we have to do
// return with energy == zero
if (!is_enabled(f_cv_active)) return 0.;
// add the biases' force, which at this point should already have
// been summed over each bias using this colvar
@ -1236,7 +1434,24 @@ cvm::real colvar::update_forces_energy()
cvm::log("Updating extended-Lagrangian degree of freedom.\n");
}
cvm::real dt = cvm::dt();
if (prev_timestep > -1) {
// Keep track of slow timestep to integrate MTS colvars
// the colvar checks the interval after waking up twice
int n_timesteps = cvm::step_relative() - prev_timestep;
if (n_timesteps != 0 && n_timesteps != time_step_factor) {
cvm::error("Error: extended-Lagrangian " + description + " has timeStepFactor " +
cvm::to_str(time_step_factor) + ", but was activated after " + cvm::to_str(n_timesteps) +
" steps at timestep " + cvm::to_str(cvm::step_absolute()) + " (relative step: " +
cvm::to_str(cvm::step_relative()) + ").\n" +
"Make sure that this colvar is requested by biases at multiples of timeStepFactor.\n");
return 0.;
}
}
prev_timestep = cvm::step_relative();
// Integrate with slow timestep (if time_step_factor != 1)
cvm::real dt = cvm::dt() * cvm::real(time_step_factor);
colvarvalue f_ext(fr.type()); // force acting on the extended variable
f_ext.reset();
@ -1248,18 +1463,17 @@ cvm::real colvar::update_forces_energy()
// - after this code block, colvar force to be applied to atomic coordinates
// ie. spring force (fb_actual will be added just below)
fr = f;
f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
// External force has been scaled for a 1-timestep impulse, scale it back because we will
// integrate it with the colvar's own timestep factor
f_ext = f / cvm::real(time_step_factor);
f_ext += (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
f = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);
// Coupling force is a slow force, to be applied to atomic coords impulse-style
f *= cvm::real(time_step_factor);
if (is_enabled(f_cv_subtract_applied_force)) {
// Report a "system" force without the biases on this colvar
// that is, just the spring force
ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
} else {
// The total force acting on the extended variable is f_ext
// This will be used in the next timestep
ft_reported = f_ext;
}
// leapfrog: starting from x_i, f_i, v_(i-1/2)
vr += (0.5 * dt) * f_ext / ext_mass;
@ -1279,13 +1493,10 @@ cvm::real colvar::update_forces_energy()
if (this->is_enabled(f_cv_periodic)) this->wrap(xr);
}
// Now adding the force on the actual colvar (for those biases who
// Now adding the force on the actual colvar (for those biases that
// bypass the extended Lagrangian mass)
f += fb_actual;
// Store force to be applied, possibly summed over several timesteps
f_accumulated += f;
if (is_enabled(f_cv_fdiff_velocity)) {
// set it for the next step
x_old = x;
@ -1306,7 +1517,7 @@ void colvar::communicate_forces()
size_t i;
if (cvm::debug()) {
cvm::log("Communicating forces from colvar \""+this->name+"\".\n");
cvm::log("Force to be applied: " + cvm::to_str(f_accumulated) + "\n");
cvm::log("Force to be applied: " + cvm::to_str(f) + "\n");
}
if (is_enabled(f_cv_scripted)) {
@ -1333,14 +1544,42 @@ void colvar::communicate_forces()
if (!cvcs[i]->is_enabled()) continue;
// cvc force is colvar force times colvar/cvc Jacobian
// (vector-matrix product)
(cvcs[i])->apply_force(colvarvalue(f_accumulated.as_vector() * func_grads[grad_index++],
(cvcs[i])->apply_force(colvarvalue(f.as_vector() * func_grads[grad_index++],
cvcs[i]->value().type()));
}
#ifdef LEPTON
} else if (is_enabled(f_cv_custom_function)) {
size_t r = 0; // index in the vector of variable references
size_t e = 0; // index of the gradient evaluator
for (size_t i = 0; i < cvcs.size(); i++) { // gradient with respect to cvc i
cvm::matrix2d<cvm::real> jacobian (x.size(), cvcs[i]->value().size());
for (size_t j = 0; j < cvcs[i]->value().size(); j++) { // j-th element
for (size_t c = 0; c < x.size(); c++) { // derivative of scalar element c of the colvarvalue
// Feed cvc values to the evaluator
for (size_t k = 0; k < cvcs.size(); k++) { //
for (size_t l = 0; l < cvcs[k]->value().size(); l++) {
*(grad_eval_var_refs[r++]) = cvcs[k]->value()[l];
}
}
jacobian[c][j] = gradient_evaluators[e++]->evaluate();
}
}
// cvc force is colvar force times colvar/cvc Jacobian
// (vector-matrix product)
(cvcs[i])->apply_force(colvarvalue(f.as_vector() * jacobian,
cvcs[i]->value().type()));
}
#endif
} else if (x.type() == colvarvalue::type_scalar) {
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->is_enabled()) continue;
(cvcs[i])->apply_force(f_accumulated * (cvcs[i])->sup_coeff *
(cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff *
cvm::real((cvcs[i])->sup_np) *
(std::pow((cvcs[i])->value().real_value,
(cvcs[i])->sup_np-1)) );
@ -1350,14 +1589,10 @@ void colvar::communicate_forces()
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->is_enabled()) continue;
(cvcs[i])->apply_force(f_accumulated * (cvcs[i])->sup_coeff);
(cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff);
}
}
// Accumulated forces have been applied, impulse-style
// Reset to start accumulating again
f_accumulated.reset();
if (cvm::debug())
cvm::log("Done communicating forces from colvar \""+this->name+"\".\n");
}
@ -1394,7 +1629,7 @@ int colvar::update_cvc_flags()
cvm::error("ERROR: All CVCs are disabled for colvar " + this->name +"\n");
return COLVARS_ERROR;
}
cvc_flags.resize(0);
cvc_flags.clear();
}
return COLVARS_OK;
@ -1744,16 +1979,15 @@ int colvar::write_output_files()
cvm::log("Writing acf to file \""+acf_outfile+"\".\n");
cvm::backup_file(acf_outfile.c_str());
cvm::ofstream acf_os(acf_outfile.c_str());
if (! acf_os.is_open()) {
cvm::error("Cannot open file \""+acf_outfile+"\".\n", FILE_ERROR);
}
write_acf(acf_os);
acf_os.close();
std::ostream *acf_os = cvm::proxy->output_stream(acf_outfile);
if (!acf_os) return cvm::get_error();
write_acf(*acf_os);
cvm::proxy->close_output_stream(acf_outfile);
}
if (runave_os.is_open()) {
runave_os.close();
if (runave_os) {
cvm::proxy->close_output_stream(runave_outfile);
runave_os = NULL;
}
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
@ -2031,7 +2265,7 @@ void colvar::calc_runave()
}
runave_variance *= 1.0 / cvm::real(runave_length-1);
runave_os << std::setw(cvm::it_width) << cvm::step_relative()
*runave_os << std::setw(cvm::it_width) << cvm::step_relative()
<< " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< runave << " "

47
lib/colvars/colvar.h
View File

@ -19,6 +19,9 @@
#include "colvarparse.h"
#include "colvardeps.h"
#ifdef LEPTON
#include "Lepton.h" // for runtime custom expressions
#endif
/// \brief A collective variable (main class); to be defined, it needs
/// at least one object of a derived class of colvar::cvc; it
@ -89,7 +92,10 @@ public:
return cv_features;
}
int refresh_deps();
/// Implements possible actions to be carried out
/// when a given feature is enabled
/// This overloads the base function in colvardeps
void do_feature_side_effects(int id);
/// List of biases that depend on this colvar
std::vector<colvarbias *> biases;
@ -235,6 +241,9 @@ public:
/// Parse the CVC configuration and allocate their data
int init_components(std::string const &conf);
/// Parse parameters for custom function with Lepton
int init_custom_function(std::string const &conf);
/// Init defaults for grid options
int init_grid_parameters(std::string const &conf);
@ -334,24 +343,13 @@ protected:
/// Sum of square coefficients for active cvcs
cvm::real active_cvc_square_norm;
/// Time step multiplier (for coarse-time-step colvars)
/// Colvar will only be calculated at those times; biases may ignore the information and
/// always update their own forces (which is typically inexpensive) especially if
/// they rely on other colvars. In this case, the colvar will accumulate forces applied between
/// colvar updates. Alternately they may use it to calculate "impulse" biasing
/// forces at longer intervals. Impulse forces must be multiplied by the timestep factor.
int time_step_factor;
/// Biasing force collected between updates, to be applied at next update for coarse-time-step colvars
colvarvalue f_accumulated;
/// \brief Absolute timestep number when this colvar was last updated
int prev_timestep;
public:
/// \brief Return the number of CVC objects with an active flag (as set by update_cvc_flags)
inline size_t num_active_cvcs() const { return n_active_cvcs; }
/// \brief returns time_step_factor
inline int get_time_step_factor() const {return time_step_factor;}
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to calculate square distances and gradients
///
@ -484,7 +482,9 @@ protected:
/// Timesteps to skip between two values in the running average series
size_t runave_stride;
/// Name of the file to write the running average
cvm::ofstream runave_os;
std::string runave_outfile;
/// File to write the running average
std::ostream *runave_os;
/// Current value of the running average
colvarvalue runave;
/// Current value of the square deviation from the running average
@ -508,6 +508,8 @@ public:
class distance;
class distance_z;
class distance_xy;
class polar_theta;
class polar_phi;
class distance_inv;
class distance_pairs;
class angle;
@ -556,6 +558,21 @@ private:
/// when using scriptedFunction
std::vector<const colvarvalue *> sorted_cvc_values;
#ifdef LEPTON
/// Vector of evaluators for custom functions using Lepton
std::vector<Lepton::CompiledExpression *> value_evaluators;
/// Vector of evaluators for gradients of custom functions
std::vector<Lepton::CompiledExpression *> gradient_evaluators;
/// Vector of references to cvc values to be passed to Lepton evaluators
std::vector<double *> value_eval_var_refs;
std::vector<double *> grad_eval_var_refs;
/// Unused value that is written to when a variable simplifies out of a Lepton expression
double dev_null;
#endif
public:
/// \brief Sorted array of (zero-based) IDs for all atoms involved
std::vector<int> atom_ids;

214
lib/colvars/colvaratoms.cpp
View File

@ -67,18 +67,16 @@ cvm::atom::~atom()
// TODO change this arrangement
// Note: "conf" is the configuration of the cvc who is using this atom group;
// "key" is the name of the atom group (e.g. "atoms", "group1", "group2", ...)
cvm::atom_group::atom_group(std::string const &conf,
char const *key_in)
cvm::atom_group::atom_group()
{
init();
}
cvm::atom_group::atom_group(char const *key_in)
{
key = key_in;
cvm::log("Defining atom group \"" + key + "\".\n");
init();
// real work is done by parse
parse(conf);
setup();
}
@ -90,12 +88,6 @@ cvm::atom_group::atom_group(std::vector<cvm::atom> const &atoms_in)
}
cvm::atom_group::atom_group()
{
init();
}
cvm::atom_group::~atom_group()
{
if (is_enabled(f_ag_scalable) && !b_dummy) {
@ -180,7 +172,7 @@ int cvm::atom_group::init()
{
if (!key.size()) key = "unnamed";
description = "atom group " + key;
// These will be overwritten by parse(), if initializing from a config string
// These may be overwritten by parse(), if a name is provided
atoms.clear();
@ -193,7 +185,6 @@ int cvm::atom_group::init()
b_center = false;
b_rotate = false;
b_user_defined_fit = false;
b_fit_gradients = false;
fitting_group = NULL;
noforce = false;
@ -265,34 +256,10 @@ void cvm::atom_group::update_total_charge()
}
int cvm::atom_group::parse(std::string const &conf)
int cvm::atom_group::parse(std::string const &group_conf)
{
std::string group_conf;
// TODO move this to the cvc class constructor/init
// save_delimiters is set to false for this call, because "conf" is
// not the config string of this group, but of its parent object
// (which has already taken care of the delimiters)
save_delimiters = false;
key_lookup(conf, key.c_str(), group_conf, dummy_pos);
// restoring the normal value, because we do want keywords checked
// inside "group_conf"
save_delimiters = true;
if (group_conf.size() == 0) {
cvm::error("Error: atom group \""+key+
"\" is set, but has no definition.\n",
INPUT_ERROR);
return COLVARS_ERROR;
}
cvm::increase_depth();
cvm::log("Initializing atom group \""+key+"\".\n");
description = "atom group " + key;
// whether or not to include messages in the log
// colvarparse::Parse_Mode mode = parse_silent;
// {
@ -304,10 +271,53 @@ int cvm::atom_group::parse(std::string const &conf)
int parse_error = COLVARS_OK;
// Optional group name will let other groups reuse atom definition
if (get_keyval(group_conf, "name", name)) {
if ((cvm::atom_group_by_name(this->name) != NULL) &&
(cvm::atom_group_by_name(this->name) != this)) {
cvm::error("Error: this atom group cannot have the same name, \""+this->name+
"\", as another atom group.\n",
INPUT_ERROR);
return INPUT_ERROR;
}
cvm::main()->register_named_atom_group(this);
description = "atom group " + name;
}
// We need to know about fitting to decide whether the group is scalable
// and we need to know about scalability before adding atoms
bool b_defined_center = get_keyval(group_conf, "centerReference", b_center, false);
bool b_defined_rotate = get_keyval(group_conf, "rotateReference", b_rotate, false);
// is the user setting explicit options?
b_user_defined_fit = b_defined_center || b_defined_rotate;
if (is_available(f_ag_scalable_com) && !b_rotate && !b_center) {
enable(f_ag_scalable_com);
enable(f_ag_scalable);
}
{
std::string atoms_of = "";
if (get_keyval(group_conf, "atomsOfGroup", atoms_of)) {
atom_group * ag = atom_group_by_name(atoms_of);
if (ag == NULL) {
cvm::error("Error: cannot find atom group with name " + atoms_of + ".\n");
return COLVARS_ERROR;
}
parse_error |= add_atoms_of_group(ag);
}
}
// if (get_keyval(group_conf, "copyOfGroup", source)) {
// // Goal: Initialize this as a full copy
// // for this we'll need an atom_group copy constructor
// return COLVARS_OK;
// }
{
std::string numbers_conf = "";
size_t pos = 0;
while (key_lookup(group_conf, "atomNumbers", numbers_conf, pos)) {
while (key_lookup(group_conf, "atomNumbers", &numbers_conf, &pos)) {
parse_error |= add_atom_numbers(numbers_conf);
numbers_conf = "";
}
@ -325,7 +335,7 @@ int cvm::atom_group::parse(std::string const &conf)
std::string range_conf = "";
size_t pos = 0;
while (key_lookup(group_conf, "atomNumbersRange",
range_conf, pos)) {
&range_conf, &pos)) {
parse_error |= add_atom_numbers_range(range_conf);
range_conf = "";
}
@ -347,7 +357,7 @@ int cvm::atom_group::parse(std::string const &conf)
size_t range_count = 0;
psii = psf_segids.begin();
while (key_lookup(group_conf, "atomNameResidueRange",
range_conf, pos)) {
&range_conf, &pos)) {
range_count++;
if (psf_segids.size() && (range_count > psf_segids.size())) {
cvm::error("Error: more instances of \"atomNameResidueRange\" than "
@ -415,14 +425,9 @@ int cvm::atom_group::parse(std::string const &conf)
}
}
// We need to know the fitting options to decide whether the group is scalable
// Now that atoms are defined we can parse the detailed fitting options
parse_error |= parse_fitting_options(group_conf);
if (is_available(f_ag_scalable_com) && !b_rotate && !b_center) {
enable(f_ag_scalable_com);
enable(f_ag_scalable);
}
if (is_enabled(f_ag_scalable) && !b_dummy) {
cvm::log("Enabling scalable calculation for group \""+this->key+"\".\n");
index = (cvm::proxy)->init_atom_group(atoms_ids);
@ -431,13 +436,6 @@ int cvm::atom_group::parse(std::string const &conf)
bool b_print_atom_ids = false;
get_keyval(group_conf, "printAtomIDs", b_print_atom_ids, false, colvarparse::parse_silent);
// TODO move this to colvarparse object
check_keywords(group_conf, key.c_str());
if (cvm::get_error()) {
cvm::error("Error setting up atom group \""+key+"\".");
return COLVARS_ERROR;
}
// Calculate all required properties (such as total mass)
setup();
@ -454,12 +452,41 @@ int cvm::atom_group::parse(std::string const &conf)
cvm::log(print_atom_ids());
}
cvm::decrease_depth();
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int cvm::atom_group::add_atoms_of_group(atom_group const * ag)
{
std::vector<int> const &source_ids = ag->atoms_ids;
if (source_ids.size()) {
atoms_ids.reserve(atoms_ids.size()+source_ids.size());
if (is_enabled(f_ag_scalable)) {
for (size_t i = 0; i < source_ids.size(); i++) {
add_atom_id(source_ids[i]);
}
} else {
atoms.reserve(atoms.size()+source_ids.size());
for (size_t i = 0; i < source_ids.size(); i++) {
// We could use the atom copy constructor, but only if the source
// group is not scalable - whereas this works in both cases
// atom constructor expects 1-based atom number
add_atom(cvm::atom(source_ids[i] + 1));
}
}
if (cvm::get_error()) return COLVARS_ERROR;
} else {
cvm::error("Error: source atom group contains no atoms\".\n", INPUT_ERROR);
return COLVARS_ERROR;
}
return COLVARS_OK;
}
int cvm::atom_group::add_atom_numbers(std::string const &numbers_conf)
{
std::vector<int> atom_indexes;
@ -629,13 +656,6 @@ std::string const cvm::atom_group::print_atom_ids() const
int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
{
bool b_defined_center = get_keyval(group_conf, "centerReference", b_center, false);
bool b_defined_rotate = get_keyval(group_conf, "rotateReference", b_rotate, false);
// is the user setting explicit options?
b_user_defined_fit = b_defined_center || b_defined_rotate;
get_keyval(group_conf, "enableFitGradients", b_fit_gradients, true);
if (b_center || b_rotate) {
if (b_dummy)
@ -643,27 +663,31 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
"cannot be defined for a dummy atom.\n");
bool b_ref_pos_group = false;
if (key_lookup(group_conf, "refPositionsGroup")) {
std::string fitting_group_conf;
if (key_lookup(group_conf, "refPositionsGroup", &fitting_group_conf)) {
b_ref_pos_group = true;
cvm::log("Warning: keyword \"refPositionsGroup\" is deprecated: please use \"fittingGroup\" instead.\n");
}
if (b_ref_pos_group || key_lookup(group_conf, "fittingGroup")) {
if (b_ref_pos_group || key_lookup(group_conf, "fittingGroup", &fitting_group_conf)) {
// instead of this group, define another group to compute the fit
if (fitting_group) {
cvm::error("Error: the atom group \""+
key+"\" has already a reference group "
"for the rototranslational fit, which was communicated by the "
"colvar component. You should not use fittingGroup "
"in this case.\n");
"in this case.\n", INPUT_ERROR);
return INPUT_ERROR;
}
cvm::log("Within atom group \""+key+"\":\n");
fitting_group = b_ref_pos_group ?
new atom_group(group_conf, "refPositionsGroup") :
new atom_group(group_conf, "fittingGroup");
// regardless of the configuration, fit gradients must be calculated by fittingGroup
fitting_group->b_fit_gradients = this->b_fit_gradients;
fitting_group = new atom_group("fittingGroup");
if (fitting_group->parse(fitting_group_conf) == COLVARS_OK) {
fitting_group->check_keywords(fitting_group_conf, "fittingGroup");
if (cvm::get_error()) {
cvm::error("Error setting up atom group \"fittingGroup\".", INPUT_ERROR);
return INPUT_ERROR;
}
}
}
atom_group *group_for_fit = fitting_group ? fitting_group : this;
@ -720,11 +744,6 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
return COLVARS_ERROR;
}
if (b_fit_gradients) {
group_for_fit->fit_gradients.assign(group_for_fit->size(), cvm::atom_pos(0.0, 0.0, 0.0));
rot.request_group1_gradients(group_for_fit->size());
}
if (b_rotate && !noforce) {
cvm::log("Warning: atom group \""+key+
"\" will be aligned to a fixed orientation given by the reference positions provided. "
@ -737,10 +756,37 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
}
}
// Enable fit gradient calculation only if necessary, and not disabled by the user
// This must happen after fitting group is defined so that side-effects are performed
// properly (ie. allocating fitting group gradients)
{
bool b_fit_gradients;
get_keyval(group_conf, "enableFitGradients", b_fit_gradients, true);
if (b_fit_gradients && (b_center || b_rotate)) {
enable(f_ag_fit_gradients);
}
}
return COLVARS_OK;
}
void cvm::atom_group::do_feature_side_effects(int id)
{
// If enabled features are changed upstream, the features below should be refreshed
switch (id) {
case f_ag_fit_gradients:
if (b_center || b_rotate) {
atom_group *group_for_fit = fitting_group ? fitting_group : this;
group_for_fit->fit_gradients.assign(group_for_fit->size(), cvm::atom_pos(0.0, 0.0, 0.0));
rot.request_group1_gradients(group_for_fit->size());
}
break;
}
}
int cvm::atom_group::create_sorted_ids(void)
{
// Only do the work if the vector is not yet populated
@ -1000,12 +1046,12 @@ void cvm::atom_group::set_weighted_gradient(cvm::rvector const &grad)
void cvm::atom_group::calc_fit_gradients()
{
if (b_dummy) return;
if (b_dummy || ! is_enabled(f_ag_fit_gradients)) return;
if (cvm::debug())
cvm::log("Calculating fit gradients.\n");
atom_group *group_for_fit = fitting_group ? fitting_group : this;
cvm::atom_group *group_for_fit = fitting_group ? fitting_group : this;
if (b_center) {
// add the center of geometry contribution to the gradients
@ -1190,7 +1236,7 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force)
}
}
if ((b_center || b_rotate) && b_fit_gradients) {
if ((b_center || b_rotate) && is_enabled(f_ag_fit_gradients)) {
atom_group *group_for_fit = fitting_group ? fitting_group : this;

39
lib/colvars/colvaratoms.h
View File

@ -150,12 +150,21 @@ class colvarmodule::atom_group
{
public:
/// \brief Initialize the group by looking up its configuration
/// string in conf and parsing it; this is actually done by parse(),
/// which is a member function so that a group can be initialized
/// also after construction
atom_group(std::string const &conf,
char const *key);
/// \brief Default constructor
atom_group();
/// \brief Create a group object, assign a name to it
atom_group(char const *key);
/// \brief Initialize the group after a (temporary) vector of atoms
atom_group(std::vector<cvm::atom> const &atoms_in);
/// \brief Destructor
~atom_group();
/// \brief Optional name to reuse properties of this in other groups
std::string name;
/// \brief Keyword used to define the group
// TODO Make this field part of the data structures that link a group to a CVC
@ -172,15 +181,13 @@ public:
int parse(std::string const &conf);
int add_atom_numbers(std::string const &numbers_conf);
int add_atoms_of_group(atom_group const * ag);
int add_index_group(std::string const &index_group_name);
int add_atom_numbers_range(std::string const &range_conf);
int add_atom_name_residue_range(std::string const &psf_segid,
std::string const &range_conf);
int parse_fitting_options(std::string const &group_conf);
/// \brief Initialize the group after a (temporary) vector of atoms
atom_group(std::vector<cvm::atom> const &atoms_in);
/// \brief Add an atom object to this group
int add_atom(cvm::atom const &a);
@ -203,12 +210,6 @@ public:
return ag_features;
}
/// \brief Default constructor
atom_group();
/// \brief Destructor
~atom_group();
protected:
/// \brief Array of atom objects
@ -294,10 +295,6 @@ public:
/// cvc's (eg rmsd, eigenvector) will not override the user's choice
bool b_user_defined_fit;
/// \brief Whether or not the derivatives of the roto-translation
/// should be included when calculating the colvar's gradients (default: yes)
bool b_fit_gradients;
/// \brief use reference coordinates for b_center or b_rotate
std::vector<cvm::atom_pos> ref_pos;
@ -464,6 +461,10 @@ public:
/// apply_colvar_force() once that is implemented for non-scalar values
void apply_force(cvm::rvector const &force);
/// Implements possible actions to be carried out
/// when a given feature is enabled
/// This overloads the base function in colvardeps
void do_feature_side_effects(int id);
};

24
lib/colvars/colvarbias.cpp
View File

@ -23,9 +23,7 @@ colvarbias::colvarbias(char const *key)
b_output_energy = false;
reset();
state_file_step = 0;
// Start in active state by default
enable(f_cvb_active);
description = "uninitialized " + cvm::to_str(key) + " bias";
}
@ -74,7 +72,6 @@ int colvarbias::init(std::string const &conf)
cvm::error("Error: no collective variables specified.\n", INPUT_ERROR);
return INPUT_ERROR;
}
} else {
cvm::log("Reinitializing bias \""+name+"\".\n");
}
@ -83,6 +80,16 @@ int colvarbias::init(std::string const &conf)
get_keyval(conf, "outputEnergy", b_output_energy, b_output_energy);
get_keyval(conf, "timeStepFactor", time_step_factor, 1);
if (time_step_factor < 1) {
cvm::error("Error: timeStepFactor must be 1 or greater.\n");
return COLVARS_ERROR;
}
// Now that children are defined, we can solve dependencies
enable(f_cvb_active);
if (cvm::debug()) print_state();
return COLVARS_OK;
}
@ -110,6 +117,8 @@ colvarbias::~colvarbias()
int colvarbias::clear()
{
free_children_deps();
// Remove references to this bias from colvars
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
@ -200,7 +209,12 @@ void colvarbias::communicate_forces()
cvm::log("Communicating a force to colvar \""+
variables(i)->name+"\".\n");
}
variables(i)->add_bias_force(colvar_forces[i]);
// Impulse-style multiple timestep
// Note that biases with different values of time_step_factor
// may send forces to the same colvar
// which is why rescaling has to happen now: the colvar is not
// aware of this bias' time_step_factor
variables(i)->add_bias_force(cvm::real(time_step_factor) * colvar_forces[i]);
}
}

120
lib/colvars/colvarbias_abf.cpp
View File

@ -71,10 +71,17 @@ int colvarbias_abf::init(std::string const &conf)
// shared ABF
get_keyval(conf, "shared", shared_on, false);
if (shared_on) {
if (!cvm::replica_enabled() || cvm::replica_num() <= 1)
if (!cvm::replica_enabled() || cvm::replica_num() <= 1) {
cvm::error("Error: shared ABF requires more than one replica.");
else
return COLVARS_ERROR;
}
cvm::log("shared ABF will be applied among "+ cvm::to_str(cvm::replica_num()) + " replicas.\n");
if (cvm::proxy->smp_enabled() == COLVARS_OK) {
cvm::error("Error: shared ABF is currently not available with SMP parallelism; "
"please set \"SMP off\" at the top of the Colvars configuration file.\n",
COLVARS_NOT_IMPLEMENTED);
return COLVARS_NOT_IMPLEMENTED;
}
// If shared_freq is not set, we default to output_freq
get_keyval(conf, "sharedFreq", shared_freq, output_freq);
@ -84,11 +91,11 @@ int colvarbias_abf::init(std::string const &conf)
if (colvars.size() == 0) {
cvm::error("Error: no collective variables specified for the ABF bias.\n");
return COLVARS_ERROR;
}
if (update_bias) {
// Request calculation of total force (which also checks for availability)
// TODO - change this to a dependency - needs ABF-specific features
// Request calculation of total force
if(enable(f_cvb_get_total_force)) return cvm::get_error();
}
@ -108,6 +115,16 @@ int colvarbias_abf::init(std::string const &conf)
if (colvars[i]->is_enabled(f_cv_extended_Lagrangian))
b_extended = true;
// Cannot mix and match coarse time steps with ABF because it gives
// wrong total force averages - total force needs to be averaged over
// every time step
if (colvars[i]->get_time_step_factor() != time_step_factor) {
cvm::error("Error: " + colvars[i]->description + " has a value of timeStepFactor ("
+ cvm::to_str(colvars[i]->get_time_step_factor()) + ") different from that of "
+ description + " (" + cvm::to_str(time_step_factor) + ").\n");
return COLVARS_ERROR;
}
// Here we could check for orthogonality of the Cartesian coordinates
// and make it just a warning if some parameter is set?
}
@ -287,7 +304,7 @@ int colvarbias_abf::update()
// Factor that ensures smooth introduction of the force
if ( count < full_samples ) {
fact = ( count < min_samples) ? 0.0 :
fact = (count < min_samples) ? 0.0 :
(cvm::real(count - min_samples)) / (cvm::real(full_samples - min_samples));
}
@ -434,62 +451,57 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::string gradients_out_name = prefix + ".grad";
std::ios::openmode mode = (append ? std::ios::app : std::ios::out);
cvm::ofstream samples_os;
cvm::ofstream gradients_os;
if (!append) cvm::backup_file(samples_out_name.c_str());
samples_os.open(samples_out_name.c_str(), mode);
if (!samples_os.is_open()) {
std::ostream *samples_os =
cvm::proxy->output_stream(samples_out_name, mode);
if (!samples_os) {
cvm::error("Error opening ABF samples file " + samples_out_name + " for writing");
}
samples->write_multicol(samples_os);
samples_os.close();
samples->write_multicol(*samples_os);
cvm::proxy->close_output_stream(samples_out_name);
if (!append) cvm::backup_file(gradients_out_name.c_str());
gradients_os.open(gradients_out_name.c_str(), mode);
if (!gradients_os.is_open()) {
std::ostream *gradients_os =
cvm::proxy->output_stream(gradients_out_name, mode);
if (!gradients_os) {
cvm::error("Error opening ABF gradient file " + gradients_out_name + " for writing");
}
gradients->write_multicol(gradients_os);
gradients_os.close();
gradients->write_multicol(*gradients_os);
cvm::proxy->close_output_stream(gradients_out_name);
if (colvars.size() == 1) {
std::string pmf_out_name = prefix + ".pmf";
if (!append) cvm::backup_file(pmf_out_name.c_str());
cvm::ofstream pmf_os;
// Do numerical integration and output a PMF
pmf_os.open(pmf_out_name.c_str(), mode);
if (!pmf_os.is_open()) cvm::error("Error opening pmf file " + pmf_out_name + " for writing");
gradients->write_1D_integral(pmf_os);
pmf_os << std::endl;
pmf_os.close();
std::string pmf_out_name = prefix + ".pmf";
std::ostream *pmf_os = cvm::proxy->output_stream(pmf_out_name, mode);
if (!pmf_os) {
cvm::error("Error opening pmf file " + pmf_out_name + " for writing");
}
gradients->write_1D_integral(*pmf_os);
*pmf_os << std::endl;
cvm::proxy->close_output_stream(pmf_out_name);
}
if (z_gradients) {
// Write eABF-related quantities
std::string z_samples_out_name = prefix + ".zcount";
cvm::ofstream z_samples_os;
if (!append) cvm::backup_file(z_samples_out_name.c_str());
z_samples_os.open(z_samples_out_name.c_str(), mode);
if (!z_samples_os.is_open()) {
std::ostream *z_samples_os =
cvm::proxy->output_stream(z_samples_out_name, mode);
if (!z_samples_os) {
cvm::error("Error opening eABF z-histogram file " + z_samples_out_name + " for writing");
}
z_samples->write_multicol(z_samples_os);
z_samples_os.close();
z_samples->write_multicol(*z_samples_os);
cvm::proxy->close_output_stream(z_samples_out_name);
if (b_czar_window_file) {
std::string z_gradients_out_name = prefix + ".zgrad";
cvm::ofstream z_gradients_os;
if (!append) cvm::backup_file(z_gradients_out_name.c_str());
z_gradients_os.open(z_gradients_out_name.c_str(), mode);
if (!z_gradients_os.is_open()) {
std::ostream *z_gradients_os =
cvm::proxy->output_stream(z_gradients_out_name, mode);
if (!z_gradients_os) {
cvm::error("Error opening eABF z-gradient file " + z_gradients_out_name + " for writing");
}
z_gradients->write_multicol(z_gradients_os);
z_gradients_os.close();
z_gradients->write_multicol(*z_gradients_os);
cvm::proxy->close_output_stream(z_gradients_out_name);
}
// Calculate CZAR estimator of gradients
@ -503,26 +515,24 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
}
std::string czar_gradients_out_name = prefix + ".czar.grad";
cvm::ofstream czar_gradients_os;
if (!append) cvm::backup_file(czar_gradients_out_name.c_str());
czar_gradients_os.open(czar_gradients_out_name.c_str(), mode);
if (!czar_gradients_os.is_open()) {
std::ostream *czar_gradients_os =
cvm::proxy->output_stream(czar_gradients_out_name, mode);
if (!czar_gradients_os) {
cvm::error("Error opening CZAR gradient file " + czar_gradients_out_name + " for writing");
}
czar_gradients->write_multicol(czar_gradients_os);
czar_gradients_os.close();
czar_gradients->write_multicol(*czar_gradients_os);
cvm::proxy->close_output_stream(czar_gradients_out_name);
if (colvars.size() == 1) {
std::string czar_pmf_out_name = prefix + ".czar.pmf";
if (!append) cvm::backup_file(czar_pmf_out_name.c_str());
cvm::ofstream czar_pmf_os;
// Do numerical integration and output a PMF
czar_pmf_os.open(czar_pmf_out_name.c_str(), mode);
if (!czar_pmf_os.is_open()) cvm::error("Error opening CZAR pmf file " + czar_pmf_out_name + " for writing");
czar_gradients->write_1D_integral(czar_pmf_os);
czar_pmf_os << std::endl;
czar_pmf_os.close();
std::string czar_pmf_out_name = prefix + ".czar.pmf";
std::ostream *czar_pmf_os =
cvm::proxy->output_stream(czar_pmf_out_name, mode);
if (!czar_pmf_os) cvm::error("Error opening CZAR pmf file " + czar_pmf_out_name + " for writing");
czar_gradients->write_1D_integral(*czar_pmf_os);
*czar_pmf_os << std::endl;
cvm::proxy->close_output_stream(czar_pmf_out_name);
}
}
return;
@ -570,9 +580,13 @@ void colvarbias_abf::read_gradients_samples()
is.clear();
is.open(gradients_in_name.c_str());
if (!is.is_open()) cvm::error("Error opening ABF gradient file " + gradients_in_name + " for reading");
if (!is.is_open()) {
cvm::error("Error opening ABF gradient file " +
gradients_in_name + " for reading", INPUT_ERROR);
} else {
gradients->read_multicol(is, true);
is.close();
}
if (z_gradients) {
// Read eABF z-averaged data for CZAR

28
lib/colvars/colvarbias_histogram.cpp
View File

@ -86,8 +86,9 @@ int colvarbias_histogram::init(std::string const &conf)
{
std::string grid_conf;
if (key_lookup(conf, "histogramGrid", grid_conf)) {
if (key_lookup(conf, "histogramGrid", &grid_conf)) {
grid->parse_params(grid_conf);
grid->check_keywords(grid_conf, "histogramGrid");
}
}
@ -176,26 +177,27 @@ int colvarbias_histogram::write_output_files()
if (out_name.size()) {
cvm::log("Writing the histogram file \""+out_name+"\".\n");
cvm::backup_file(out_name.c_str());
cvm::ofstream grid_os(out_name.c_str());
if (!grid_os.is_open()) {
cvm::error("Error opening histogram file " + out_name + " for writing.\n", FILE_ERROR);
std::ostream *grid_os = cvm::proxy->output_stream(out_name);
if (!grid_os) {
return cvm::error("Error opening histogram file "+out_name+
" for writing.\n", FILE_ERROR);
}
// TODO add return code here
grid->write_multicol(grid_os);
grid_os.close();
grid->write_multicol(*grid_os);
cvm::proxy->close_output_stream(out_name);
}
if (out_name_dx.size()) {
cvm::log("Writing the histogram file \""+out_name_dx+"\".\n");
cvm::backup_file(out_name_dx.c_str());
cvm::ofstream grid_os(out_name_dx.c_str());
if (!grid_os.is_open()) {
cvm::error("Error opening histogram file " + out_name_dx + " for writing.\n", FILE_ERROR);
std::ostream *grid_os = cvm::proxy->output_stream(out_name_dx);
if (!grid_os) {
return cvm::error("Error opening histogram file "+out_name_dx+
" for writing.\n", FILE_ERROR);
}
// TODO add return code here
grid->write_opendx(grid_os);
grid_os.close();
grid->write_opendx(*grid_os);
cvm::proxy->close_output_stream(out_name_dx);
}
return COLVARS_OK;
}

141
lib/colvars/colvarbias_meta.cpp
View File

@ -36,6 +36,8 @@ colvarbias_meta::colvarbias_meta(char const *key)
: colvarbias(key)
{
new_hills_begin = hills.end();
hills_traj_os = NULL;
replica_hills_os = NULL;
}
@ -163,7 +165,6 @@ int colvarbias_meta::init(std::string const &conf)
cvm::log("Done initializing the metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
save_delimiters = false;
return COLVARS_OK;
}
@ -239,11 +240,15 @@ colvarbias_meta::~colvarbias_meta()
hills_energy_gradients = NULL;
}
if (replica_hills_os.is_open())
replica_hills_os.close();
if (replica_hills_os) {
cvm::proxy->close_output_stream(replica_hills_file);
replica_hills_os = NULL;
}
if (hills_traj_os.is_open())
hills_traj_os.close();
if (hills_traj_os) {
cvm::proxy->close_output_stream(hills_traj_file_name());
hills_traj_os = NULL;
}
if(target_dist) {
delete target_dist;
@ -280,9 +285,9 @@ colvarbias_meta::create_hill(colvarbias_meta::hill const &h)
}
// output to trajectory (if specified)
if (hills_traj_os.is_open()) {
hills_traj_os << (hills.back()).output_traj();
hills_traj_os.flush();
if (hills_traj_os) {
*hills_traj_os << (hills.back()).output_traj();
cvm::proxy->flush_output_stream(hills_traj_os);
}
has_data = true;
@ -312,12 +317,12 @@ colvarbias_meta::delete_hill(hill_iter &h)
}
}
if (hills_traj_os.is_open()) {
if (hills_traj_os) {
// output to the trajectory
hills_traj_os << "# DELETED this hill: "
*hills_traj_os << "# DELETED this hill: "
<< (hills.back()).output_traj()
<< "\n";
hills_traj_os.flush();
cvm::proxy->flush_output_stream(hills_traj_os);
}
return hills.erase(h);
@ -501,12 +506,12 @@ int colvarbias_meta::update_bias()
case multiple_replicas:
create_hill(hill(hill_weight*hills_scale, colvars, hill_width, replica_id));
if (replica_hills_os.is_open()) {
replica_hills_os << hills.back();
if (replica_hills_os) {
*replica_hills_os << hills.back();
} else {
cvm::fatal_error("Error: in metadynamics bias \""+this->name+"\""+
return cvm::error("Error: in metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
" while writing hills for the other replicas.\n");
" while writing hills for the other replicas.\n", FILE_ERROR);
}
break;
}
@ -904,8 +909,9 @@ int colvarbias_meta::replica_share()
// reread the replicas registry
update_replicas_registry();
// empty the output buffer
if (replica_hills_os.is_open())
replica_hills_os.flush();
if (replica_hills_os) {
cvm::proxy->flush_output_stream(replica_hills_os);
}
read_replica_files();
}
return COLVARS_OK;
@ -1421,7 +1427,7 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
// it is safer to read colvarvalue objects one at a time;
// TODO: change this it later
std::string centers_input;
key_lookup(data, "centers", centers_input);
key_lookup(data, "centers", &centers_input);
std::istringstream centers_is(centers_input);
for (size_t i = 0; i < num_variables(); i++) {
centers_is >> h_centers[i];
@ -1521,13 +1527,11 @@ int colvarbias_meta::setup_output()
// for the others to read
// open the "hills" buffer file
if (!replica_hills_os.is_open()) {
cvm::backup_file(replica_hills_file.c_str());
replica_hills_os.open(replica_hills_file.c_str());
if (!replica_hills_os.is_open())
cvm::error("Error: in opening file \""+
replica_hills_file+"\" for writing.\n", FILE_ERROR);
replica_hills_os.setf(std::ios::scientific, std::ios::floatfield);
if (!replica_hills_os) {
cvm::proxy->backup_file(replica_hills_file);
replica_hills_os = cvm::proxy->output_stream(replica_hills_file);
if (!replica_hills_os) return cvm::get_error();
replica_hills_os->setf(std::ios::scientific, std::ios::floatfield);
}
// write the state file (so that there is always one available)
@ -1539,43 +1543,49 @@ int colvarbias_meta::setup_output()
// if we're running without grids, use a growing list of "hills" files
// otherwise, just one state file and one "hills" file as buffer
std::ofstream list_os(replica_list_file.c_str(),
(use_grids ? std::ios::trunc : std::ios::app));
if (! list_os.is_open())
cvm::fatal_error("Error: in opening file \""+
replica_list_file+"\" for writing.\n");
list_os << "stateFile " << replica_state_file << "\n";
list_os << "hillsFile " << replica_hills_file << "\n";
list_os.close();
std::ostream *list_os =
cvm::proxy->output_stream(replica_list_file,
(use_grids ? std::ios_base::trunc :
std::ios_base::app));
if (!list_os) {
return cvm::get_error();
}
*list_os << "stateFile " << replica_state_file << "\n";
*list_os << "hillsFile " << replica_hills_file << "\n";
cvm::proxy->close_output_stream(replica_list_file);
// finally, if add a new record for this replica to the registry
// finally, add a new record for this replica to the registry
if (! registered_replica) {
std::ofstream reg_os(replicas_registry_file.c_str(), std::ios::app);
if (! reg_os.is_open())
cvm::error("Error: in opening file \""+
replicas_registry_file+"\" for writing.\n", FILE_ERROR);
reg_os << replica_id << " " << replica_list_file << "\n";
reg_os.close();
std::ostream *reg_os =
cvm::proxy->output_stream(replicas_registry_file,
std::ios::app);
if (!reg_os) {
return cvm::get_error();
}
*reg_os << replica_id << " " << replica_list_file << "\n";
cvm::proxy->close_output_stream(replicas_registry_file);
}
}
if (b_hills_traj) {
std::string const traj_file_name(cvm::output_prefix()+
if (!hills_traj_os) {
hills_traj_os = cvm::proxy->output_stream(hills_traj_file_name());
if (!hills_traj_os) return cvm::get_error();
}
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
std::string const colvarbias_meta::hills_traj_file_name() const
{
return std::string(cvm::output_prefix()+
".colvars."+this->name+
( (comm != single_replica) ?
("."+replica_id) :
("") )+
".hills.traj");
if (!hills_traj_os.is_open()) {
cvm::backup_file(traj_file_name.c_str());
hills_traj_os.open(traj_file_name.c_str());
}
if (!hills_traj_os.is_open())
cvm::error("Error: in opening hills output file \"" +
traj_file_name+"\".\n", FILE_ERROR);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
@ -1671,12 +1681,13 @@ void colvarbias_meta::write_pmf()
(dump_fes_save ?
"."+cvm::to_str(cvm::step_absolute()) : "") +
".pmf");
cvm::backup_file(fes_file_name.c_str());
cvm::ofstream fes_os(fes_file_name.c_str());
pmf->write_multicol(fes_os);
fes_os.close();
cvm::proxy->backup_file(fes_file_name);
std::ostream *fes_os = cvm::proxy->output_stream(fes_file_name);
pmf->write_multicol(*fes_os);
cvm::proxy->close_output_stream(fes_file_name);
}
}
if (comm != single_replica) {
// output the combined PMF from all replicas
pmf->reset();
@ -1695,10 +1706,10 @@ void colvarbias_meta::write_pmf()
(dump_fes_save ?
"."+cvm::to_str(cvm::step_absolute()) : "") +
".pmf");
cvm::backup_file(fes_file_name.c_str());
cvm::ofstream fes_os(fes_file_name.c_str());
pmf->write_multicol(fes_os);
fes_os.close();
cvm::proxy->backup_file(fes_file_name);
std::ostream *fes_os = cvm::proxy->output_stream(fes_file_name);
pmf->write_multicol(*fes_os);
cvm::proxy->close_output_stream(fes_file_name);
}
delete pmf;
@ -1769,13 +1780,11 @@ int colvarbias_meta::write_replica_state_file()
// rep_state_os.close();
// reopen the hills file
replica_hills_os.close();
cvm::backup_file(replica_hills_file.c_str());
replica_hills_os.open(replica_hills_file.c_str());
if (!replica_hills_os.is_open())
cvm::fatal_error("Error: in opening file \""+
replica_hills_file+"\" for writing.\n");
replica_hills_os.setf(std::ios::scientific, std::ios::floatfield);
cvm::proxy->close_output_stream(replica_hills_file);
cvm::proxy->backup_file(replica_hills_file);
replica_hills_os = cvm::proxy->output_stream(replica_hills_file);
if (!replica_hills_os) return cvm::get_error();
replica_hills_os->setf(std::ios::scientific, std::ios::floatfield);
return COLVARS_OK;
}

7
lib/colvars/colvarbias_meta.h
View File

@ -78,7 +78,10 @@ protected:
/// Write the hill logfile
bool b_hills_traj;
/// Logfile of hill management (creation and deletion)
cvm::ofstream hills_traj_os;
std::ostream *hills_traj_os;
/// Name of the hill logfile
std::string const hills_traj_file_name() const;
/// \brief List of hills used on this bias (total); if a grid is
/// employed, these don't need to be updated at every time step
@ -241,7 +244,7 @@ protected:
std::string replica_hills_file;
/// \brief Output stream corresponding to replica_hills_file
cvm::ofstream replica_hills_os;
std::ostream *replica_hills_os;
/// Position within replica_hills_file (when reading it)
int replica_hills_file_pos;

36
lib/colvars/colvarbias_restraint.cpp
View File

@ -853,6 +853,21 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
get_keyval(conf, "upperWallConstant", upper_wall_k,
(upper_wall_k > 0.0) ? upper_wall_k : force_k);
if (lower_wall_k * upper_wall_k > 0.0) {
for (size_t i = 0; i < num_variables(); i++) {
if (variables(i)->width != 1.0)
cvm::log("The lower and upper wall force constants for colvar \""+
variables(i)->name+
"\" will be rescaled to "+
cvm::to_str(lower_wall_k /
(variables(i)->width * variables(i)->width))+
" and "+
cvm::to_str(upper_wall_k /
(variables(i)->width * variables(i)->width))+
" according to the specified width.\n");
}
}
enable(f_cvb_scalar_variables);
size_t i;
@ -869,7 +884,7 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
if (!get_keyval(conf, "lowerWalls", lower_walls, lower_walls) &&
b_null_lower_walls) {
cvm::log("Lower walls were not provided.\n");
lower_walls.resize(0);
lower_walls.clear();
}
bool b_null_upper_walls = false;
@ -884,7 +899,7 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
if (!get_keyval(conf, "upperWalls", upper_walls, upper_walls) &&
b_null_upper_walls) {
cvm::log("Upper walls were not provided.\n");
upper_walls.resize(0);
upper_walls.clear();
}
if ((lower_walls.size() == 0) && (upper_walls.size() == 0)) {
@ -954,7 +969,8 @@ void colvarbias_restraint_harmonic_walls::communicate_forces()
cvm::log("Communicating a force to colvar \""+
variables(i)->name+"\".\n");
}
variables(i)->add_bias_force_actual_value(colvar_forces[i]);
// Impulse-style multiple timestep
variables(i)->add_bias_force_actual_value(cvm::real(time_step_factor) * colvar_forces[i]);
}
}
@ -1282,9 +1298,9 @@ int colvarbias_restraint_histogram::init(std::string const &conf)
colvarbias_restraint_histogram::~colvarbias_restraint_histogram()
{
p.resize(0);
ref_p.resize(0);
p_diff.resize(0);
p.clear();
ref_p.clear();
p_diff.clear();
}
@ -1382,14 +1398,14 @@ std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os)
{
if (b_write_histogram) {
std::string file_name(cvm::output_prefix()+"."+this->name+".hist.dat");
std::ofstream os(file_name.c_str());
os << "# " << cvm::wrap_string(variables(0)->name, cvm::cv_width)
std::ostream *os = cvm::proxy->output_stream(file_name);
*os << "# " << cvm::wrap_string(variables(0)->name, cvm::cv_width)
<< " " << "p(" << cvm::wrap_string(variables(0)->name, cvm::cv_width-3)
<< ")\n";
size_t igrid;
for (igrid = 0; igrid < p.size(); igrid++) {
cvm::real const x_grid = (lower_boundary + (igrid+1)*width);
os << " "
*os << " "
<< std::setprecision(cvm::cv_prec)
<< std::setw(cvm::cv_width)
<< x_grid
@ -1398,7 +1414,7 @@ std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os)
<< std::setw(cvm::cv_width)
<< p[igrid] << "\n";
}
os.close();
cvm::proxy->close_output_stream(file_name);
}
return os;
}

76
lib/colvars/colvarcomp.cpp
View File

@ -51,6 +51,17 @@ colvar::cvc::cvc(std::string const &conf)
get_keyval_feature((colvarparse *)this, conf, "debugGradients",
f_cvc_debug_gradient, false, parse_silent);
{
bool b_no_PBC = false;
get_keyval(conf, "forceNoPBC", b_no_PBC, false);
if (b_no_PBC) {
disable(f_cvc_pbc_minimum_image);
} else {
enable(f_cvc_pbc_minimum_image);
}
// this does not use get_keyval_feature() only for backward compatibility
}
// Attempt scalable calculations when in parallel? (By default yes, if available)
get_keyval(conf, "scalable", b_try_scalable, true);
@ -94,13 +105,15 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
bool optional)
{
cvm::atom_group *group = NULL;
std::string group_conf;
if (key_lookup(conf, group_key)) {
group = new cvm::atom_group;
group->key = group_key;
if (key_lookup(conf, group_key, &group_conf)) {
group = new cvm::atom_group(group_key);
if (b_try_scalable) {
if (is_available(f_cvc_scalable_com) && is_enabled(f_cvc_com_based)) {
if (is_available(f_cvc_scalable_com)
&& is_enabled(f_cvc_com_based)
&& !is_enabled(f_cvc_debug_gradient)) {
enable(f_cvc_scalable_com);
enable(f_cvc_scalable);
// The CVC makes the feature available;
@ -111,44 +124,51 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
// TODO check for other types of parallelism here
}
if (group->parse(conf) == COLVARS_OK) {
atom_groups.push_back(group);
} else {
cvm::error("Error parsing definition for atom group \""+
std::string(group_key)+"\".\n");
if (group_conf.size() == 0) {
cvm::error("Error: atom group \""+group->key+
"\" is set, but has no definition.\n",
INPUT_ERROR);
return group;
}
cvm::increase_depth();
if (group->parse(group_conf) == COLVARS_OK) {
register_atom_group(group);
}
group->check_keywords(group_conf, group_key);
if (cvm::get_error()) {
cvm::error("Error parsing definition for atom group \""+
std::string(group_key)+"\"\n.", INPUT_ERROR);
}
cvm::decrease_depth();
} else {
if (! optional) {
cvm::error("Error: definition for atom group \""+
std::string(group_key)+"\" not found.\n");
}
}
return group;
}
int colvar::cvc::setup()
{
size_t i;
description = "cvc " + name;
for (i = 0; i < atom_groups.size(); i++) {
add_child((colvardeps *) atom_groups[i]);
}
return COLVARS_OK;
}
colvar::cvc::~cvc()
{
free_children_deps();
remove_all_children();
for (size_t i = 0; i < atom_groups.size(); i++) {
if (atom_groups[i] != NULL) delete atom_groups[i];
}
}
void colvar::cvc::read_data()
{
size_t ig;
@ -187,14 +207,26 @@ void colvar::cvc::calc_Jacobian_derivative()
}
void colvar::cvc::debug_gradients(cvm::atom_group *group)
void colvar::cvc::calc_fit_gradients()
{
// this function should work for any scalar variable:
for (size_t ig = 0; ig < atom_groups.size(); ig++) {
atom_groups[ig]->calc_fit_gradients();
}
}
void colvar::cvc::debug_gradients()
{
// this function should work for any scalar cvc:
// the only difference will be the name of the atom group (here, "group")
// NOTE: this assumes that groups for this cvc are non-overlapping,
// since atom coordinates are modified only within the current group
if (group->b_dummy) return;
cvm::log("Debugging gradients for " + description);
for (size_t ig = 0; ig < atom_groups.size(); ig++) {
cvm::atom_group *group = atom_groups[ig];
if (group->b_dummy) continue;
cvm::rotation const rot_0 = group->rot;
cvm::rotation const rot_inv = group->rot.inverse();
@ -209,7 +241,7 @@ void colvar::cvc::debug_gradients(cvm::atom_group *group)
// print the values of the fit gradients
if (group->b_rotate || group->b_center) {
if (group->b_fit_gradients) {
if (group->is_enabled(f_ag_fit_gradients)) {
size_t j;
// fit_gradients are in the simulation frame: we should print them in the rotated frame
@ -256,7 +288,7 @@ void colvar::cvc::debug_gradients(cvm::atom_group *group)
}
}
if ((group->b_fit_gradients) && (group->fitting_group != NULL)) {
if ((group->is_enabled(f_ag_fit_gradients)) && (group->fitting_group != NULL)) {
cvm::atom_group *ref_group = group->fitting_group;
group->read_positions();
group->calc_required_properties();
@ -297,7 +329,7 @@ void colvar::cvc::debug_gradients(cvm::atom_group *group)
cvm::log("Gradient sum: " + cvm::to_str(gradient_sum) +
" Fit gradient sum: " + cvm::to_str(fit_gradient_sum) +
" Total " + cvm::to_str(gradient_sum + fit_gradient_sum));
}
return;
}

75
lib/colvars/colvarcomp.h
View File

@ -146,8 +146,11 @@ public:
/// order to apply forces
virtual void calc_gradients() = 0;
/// \brief Calculate the atomic fit gradients
void calc_fit_gradients();
/// \brief Calculate finite-difference gradients alongside the analytical ones, for each Cartesian component
virtual void debug_gradients(cvm::atom_group *group);
virtual void debug_gradients();
/// \brief Calculate the total force from the system using the
/// inverse atomic gradients
@ -228,6 +231,12 @@ public:
/// e.g. atomic gradients
std::vector<cvm::atom_group *> atom_groups;
/// \brief Store a pointer to new atom group, and list as child for dependencies
inline void register_atom_group(cvm::atom_group *ag) {
atom_groups.push_back(ag);
add_child((colvardeps *) ag);
}
/// \brief Whether or not this CVC will be computed in parallel whenever possible
bool b_try_scalable;
@ -427,15 +436,77 @@ public:
};
/// \brief Colvar component: polar coordinate phi of a group
/// (colvarvalue::type_scalar type, range [-180:180])
class colvar::polar_phi
: public colvar::cvc
{
public:
polar_phi(std::string const &conf);
polar_phi();
virtual ~polar_phi() {}
protected:
cvm::atom_group *atoms;
cvm::real r, theta, phi;
public:
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
/// Redefined to handle the 2*PI periodicity
virtual cvm::real dist2(colvarvalue const &x1,
colvarvalue const &x2) const;
/// Redefined to handle the 2*PI periodicity
virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
/// Redefined to handle the 2*PI periodicity
virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
/// Redefined to handle the 2*PI periodicity
virtual void wrap(colvarvalue &x) const;
};
/// \brief Colvar component: polar coordinate theta of a group
/// (colvarvalue::type_scalar type, range [0:180])
class colvar::polar_theta
: public colvar::cvc
{
public:
polar_theta(std::string const &conf);
polar_theta();
virtual ~polar_theta() {}
protected:
cvm::atom_group *atoms;
cvm::real r, theta, phi;
public:
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
/// Redefined to override the distance ones
virtual cvm::real dist2(colvarvalue const &x1,
colvarvalue const &x2) const;
/// Redefined to override the distance ones
virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
/// Redefined to override the distance ones
virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
};
/// \brief Colvar component: average distance between two groups of atoms, weighted as the sixth power,
/// as in NMR refinements(colvarvalue::type_scalar type, range (0:*))
class colvar::distance_inv
: public colvar::distance
: public colvar::cvc
{
protected:
/// First atom group
cvm::atom_group *group1;
/// Second atom group
cvm::atom_group *group2;
/// Components of the distance vector orthogonal to the axis
int exponent;
/// Use absolute positions, ignoring PBCs when present
bool b_no_PBC;
public:
distance_inv(std::string const &conf);
distance_inv();

189
lib/colvars/colvarcomp_angles.cpp
View File

@ -45,9 +45,9 @@ colvar::angle::angle(cvm::atom const &a1,
group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
group3 = new cvm::atom_group(std::vector<cvm::atom>(1, a3));
atom_groups.push_back(group1);
atom_groups.push_back(group2);
atom_groups.push_back(group3);
register_atom_group(group1);
register_atom_group(group2);
register_atom_group(group3);
x.type(colvarvalue::type_scalar);
}
@ -66,9 +66,13 @@ void colvar::angle::calc_value()
cvm::atom_pos const g2_pos = group2->center_of_mass();
cvm::atom_pos const g3_pos = group3->center_of_mass();
r21 = cvm::position_distance(g2_pos, g1_pos);
r21 = is_enabled(f_cvc_pbc_minimum_image) ?
cvm::position_distance(g2_pos, g1_pos) :
g1_pos - g2_pos;
r21l = r21.norm();
r23 = cvm::position_distance(g2_pos, g3_pos);
r23 = is_enabled(f_cvc_pbc_minimum_image) ?
cvm::position_distance(g2_pos, g3_pos) :
g3_pos - g2_pos;
r23l = r23.norm();
cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
@ -166,9 +170,9 @@ colvar::dipole_angle::dipole_angle(cvm::atom const &a1,
group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
group3 = new cvm::atom_group(std::vector<cvm::atom>(1, a3));
atom_groups.push_back(group1);
atom_groups.push_back(group2);
atom_groups.push_back(group3);
register_atom_group(group1);
register_atom_group(group2);
register_atom_group(group3);
x.type(colvarvalue::type_scalar);
}
@ -191,7 +195,9 @@ void colvar::dipole_angle::calc_value()
r21 = group1->dipole();
r21l = r21.norm();
r23 = cvm::position_distance(g2_pos, g3_pos);
r23 = is_enabled(f_cvc_pbc_minimum_image) ?
cvm::position_distance(g2_pos, g3_pos) :
g3_pos - g2_pos;
r23l = r23.norm();
cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
@ -293,10 +299,10 @@ colvar::dihedral::dihedral(cvm::atom const &a1,
group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
group3 = new cvm::atom_group(std::vector<cvm::atom>(1, a3));
group4 = new cvm::atom_group(std::vector<cvm::atom>(1, a4));
atom_groups.push_back(group1);
atom_groups.push_back(group2);
atom_groups.push_back(group3);
atom_groups.push_back(group4);
register_atom_group(group1);
register_atom_group(group2);
register_atom_group(group3);
register_atom_group(group4);
x.type(colvarvalue::type_scalar);
@ -324,9 +330,15 @@ void colvar::dihedral::calc_value()
cvm::atom_pos const g4_pos = group4->center_of_mass();
// Usual sign convention: r12 = r2 - r1
r12 = cvm::position_distance(g1_pos, g2_pos);
r23 = cvm::position_distance(g2_pos, g3_pos);
r34 = cvm::position_distance(g3_pos, g4_pos);
r12 = is_enabled(f_cvc_pbc_minimum_image) ?
cvm::position_distance(g1_pos, g2_pos) :
g2_pos - g1_pos;
r23 = is_enabled(f_cvc_pbc_minimum_image) ?
cvm::position_distance(g2_pos, g3_pos) :
g3_pos - g2_pos;
r34 = is_enabled(f_cvc_pbc_minimum_image) ?
cvm::position_distance(g3_pos, g4_pos) :
g4_pos - g3_pos;
cvm::rvector const n1 = cvm::rvector::outer(r12, r23);
cvm::rvector const n2 = cvm::rvector::outer(r23, r34);
@ -510,3 +522,148 @@ void colvar::dihedral::wrap(colvarvalue &x) const
return;
}
colvar::polar_theta::polar_theta(std::string const &conf)
: cvc(conf)
{
function_type = "polar_theta";
enable(f_cvc_com_based);
atoms = parse_group(conf, "atoms");
init_total_force_params(conf);
x.type(colvarvalue::type_scalar);
}
colvar::polar_theta::polar_theta()
{
function_type = "polar_theta";
x.type(colvarvalue::type_scalar);
}
void colvar::polar_theta::calc_value()
{
cvm::rvector pos = atoms->center_of_mass();
r = atoms->center_of_mass().norm();
// Internal values of theta and phi are radians
theta = (r > 0.) ? std::acos(pos.z / r) : 0.;
phi = std::atan2(pos.y, pos.x);
x.real_value = (180.0/PI) * theta;
}
void colvar::polar_theta::calc_gradients()
{
if (r == 0.)
atoms->set_weighted_gradient(cvm::rvector(0., 0., 0.));
else
atoms->set_weighted_gradient(cvm::rvector(
(180.0/PI) * std::cos(theta) * std::cos(phi) / r,
(180.0/PI) * std::cos(theta) * std::sin(phi) / r,
(180.0/PI) * -std::sin(theta) / r));
}
void colvar::polar_theta::apply_force(colvarvalue const &force)
{
if (!atoms->noforce)
atoms->apply_colvar_force(force.real_value);
}
simple_scalar_dist_functions(polar_theta)
colvar::polar_phi::polar_phi(std::string const &conf)
: cvc(conf)
{
function_type = "polar_phi";
period = 360.0;
enable(f_cvc_com_based);
atoms = parse_group(conf, "atoms");
init_total_force_params(conf);
x.type(colvarvalue::type_scalar);
}
colvar::polar_phi::polar_phi()
{
function_type = "polar_phi";
period = 360.0;
x.type(colvarvalue::type_scalar);
}
void colvar::polar_phi::calc_value()
{
cvm::rvector pos = atoms->center_of_mass();
r = atoms->center_of_mass().norm();
// Internal values of theta and phi are radians
theta = (r > 0.) ? std::acos(pos.z / r) : 0.;
phi = std::atan2(pos.y, pos.x);
x.real_value = (180.0/PI) * phi;
}
void colvar::polar_phi::calc_gradients()
{
atoms->set_weighted_gradient(cvm::rvector(
(180.0/PI) * -std::sin(phi) / (r*std::sin(theta)),
(180.0/PI) * std::cos(phi) / (r*std::sin(theta)),
0.));
}
void colvar::polar_phi::apply_force(colvarvalue const &force)
{
if (!atoms->noforce)
atoms->apply_colvar_force(force.real_value);
}
// Same as dihedral, for polar_phi
cvm::real colvar::polar_phi::dist2(colvarvalue const &x1,
colvarvalue const &x2) const
{
cvm::real diff = x1.real_value - x2.real_value;
diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
return diff * diff;
}
colvarvalue colvar::polar_phi::dist2_lgrad(colvarvalue const &x1,
colvarvalue const &x2) const
{
cvm::real diff = x1.real_value - x2.real_value;
diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
return 2.0 * diff;
}
colvarvalue colvar::polar_phi::dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const
{
cvm::real diff = x1.real_value - x2.real_value;
diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
return (-2.0) * diff;
}
void colvar::polar_phi::wrap(colvarvalue &x) const
{
if ((x.real_value - wrap_center) >= 180.0) {
x.real_value -= 360.0;
return;
}
if ((x.real_value - wrap_center) < -180.0) {
x.real_value += 360.0;
return;
}
return;
}

47
lib/colvars/colvarcomp_coordnums.cpp
View File

@ -87,8 +87,10 @@ colvar::coordnum::coordnum(std::string const &conf)
group1 = parse_group(conf, "group1");
group2 = parse_group(conf, "group2");
if (group1->b_dummy)
cvm::fatal_error("Error: only group2 is allowed to be a dummy atom\n");
if (group1->b_dummy) {
cvm::error("Error: only group2 is allowed to be a dummy atom\n");
return;
}
bool const b_isotropic = get_keyval(conf, "cutoff", r0,
cvm::real(4.0 * cvm::unit_angstrom()));
@ -99,6 +101,7 @@ colvar::coordnum::coordnum(std::string const &conf)
if (b_isotropic) {
cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" at the same time.\n",
INPUT_ERROR);
return;
}
b_anisotropic = true;
@ -115,6 +118,10 @@ colvar::coordnum::coordnum(std::string const &conf)
cvm::error("Error: odd exponents provided, can only use even ones.\n", INPUT_ERROR);
}
if (!is_enabled(f_cvc_pbc_minimum_image)) {
cvm::log("Warning: only minimum-image distances are used by this variable.\n");
}
get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2->b_dummy);
}
@ -228,12 +235,13 @@ colvar::h_bond::h_bond(std::string const &conf)
get_keyval(conf, "donor", d_num, -1);
if ( (a_num == -1) || (d_num == -1) ) {
cvm::fatal_error("Error: either acceptor or donor undefined.\n");
cvm::error("Error: either acceptor or donor undefined.\n");
return;
}
cvm::atom acceptor = cvm::atom(a_num);
cvm::atom donor = cvm::atom(d_num);
atom_groups.push_back(new cvm::atom_group);
register_atom_group(new cvm::atom_group);
atom_groups[0]->add_atom(acceptor);
atom_groups[0]->add_atom(donor);
@ -242,7 +250,8 @@ colvar::h_bond::h_bond(std::string const &conf)
get_keyval(conf, "expDenom", ed, 8);
if ( (en%2) || (ed%2) ) {
cvm::fatal_error("Error: odd exponents provided, can only use even ones.\n");
cvm::error("Error: odd exponents provided, can only use even ones.\n");
return;
}
if (cvm::debug())
@ -258,7 +267,7 @@ colvar::h_bond::h_bond(cvm::atom const &acceptor,
function_type = "h_bond";
x.type(colvarvalue::type_scalar);
atom_groups.push_back(new cvm::atom_group);
register_atom_group(new cvm::atom_group);
atom_groups[0]->add_atom(acceptor);
atom_groups[0]->add_atom(donor);
}
@ -313,7 +322,12 @@ colvar::selfcoordnum::selfcoordnum(std::string const &conf)
get_keyval(conf, "expDenom", ed, int(12));
if ( (en%2) || (ed%2) ) {
cvm::fatal_error("Error: odd exponents provided, can only use even ones.\n");
cvm::error("Error: odd exponents provided, can only use even ones.\n");
return;
}
if (!is_enabled(f_cvc_pbc_minimum_image)) {
cvm::log("Warning: only minimum-image distances are used by this variable.\n");
}
}
@ -364,8 +378,10 @@ colvar::groupcoordnum::groupcoordnum(std::string const &conf)
x.type(colvarvalue::type_scalar);
// group1 and group2 are already initialized by distance()
if (group1->b_dummy || group2->b_dummy)
cvm::fatal_error("Error: neither group can be a dummy atom\n");
if (group1->b_dummy || group2->b_dummy) {
cvm::error("Error: neither group can be a dummy atom\n");
return;
}
bool const b_scale = get_keyval(conf, "cutoff", r0,
cvm::real(4.0 * cvm::unit_angstrom()));
@ -373,9 +389,11 @@ colvar::groupcoordnum::groupcoordnum(std::string const &conf)
if (get_keyval(conf, "cutoff3", r0_vec,
cvm::rvector(4.0, 4.0, 4.0), parse_silent)) {
if (b_scale)
cvm::fatal_error("Error: cannot specify \"scale\" and "
if (b_scale) {
cvm::error("Error: cannot specify \"scale\" and "
"\"scale3\" at the same time.\n");
return;
}
b_anisotropic = true;
// remove meaningless negative signs
if (r0_vec.x < 0.0) r0_vec.x *= -1.0;
@ -387,7 +405,12 @@ colvar::groupcoordnum::groupcoordnum(std::string const &conf)
get_keyval(conf, "expDenom", ed, int(12));
if ( (en%2) || (ed%2) ) {
cvm::fatal_error("Error: odd exponents provided, can only use even ones.\n");
cvm::error("Error: odd exponents provided, can only use even ones.\n");
return;
}
if (!is_enabled(f_cvc_pbc_minimum_image)) {
cvm::log("Warning: only minimum-image distances are used by this variable.\n");
}
}

91
lib/colvars/colvarcomp_distances.cpp
View File

@ -28,10 +28,6 @@ colvar::distance::distance(std::string const &conf)
group1 = parse_group(conf, "group1");
group2 = parse_group(conf, "group2");
if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) {
cvm::log("Computing distance using absolute positions (not minimal-image)");
}
init_total_force_params(conf);
x.type(colvarvalue::type_scalar);
@ -45,14 +41,13 @@ colvar::distance::distance()
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
enable(f_cvc_com_based);
b_no_PBC = false;
x.type(colvarvalue::type_scalar);
}
void colvar::distance::calc_value()
{
if (b_no_PBC) {
if (!is_enabled(f_cvc_pbc_minimum_image)) {
dist_v = group2->center_of_mass() - group1->center_of_mass();
} else {
dist_v = cvm::position_distance(group1->center_of_mass(),
@ -107,6 +102,7 @@ colvar::distance_vec::distance_vec(std::string const &conf)
{
function_type = "distance_vec";
enable(f_cvc_com_based);
enable(f_cvc_implicit_gradient);
x.type(colvarvalue::type_3vector);
}
@ -116,13 +112,14 @@ colvar::distance_vec::distance_vec()
{
function_type = "distance_vec";
enable(f_cvc_com_based);
enable(f_cvc_implicit_gradient);
x.type(colvarvalue::type_3vector);
}
void colvar::distance_vec::calc_value()
{
if (b_no_PBC) {
if (!is_enabled(f_cvc_pbc_minimum_image)) {
x.rvector_value = group2->center_of_mass() - group1->center_of_mass();
} else {
x.rvector_value = cvm::position_distance(group1->center_of_mass(),
@ -214,10 +211,6 @@ colvar::distance_z::distance_z(std::string const &conf)
fixed_axis = true;
}
if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) {
cvm::log("Computing distance using absolute positions (not minimal-image)");
}
init_total_force_params(conf);
}
@ -236,7 +229,7 @@ colvar::distance_z::distance_z()
void colvar::distance_z::calc_value()
{
if (fixed_axis) {
if (b_no_PBC) {
if (!is_enabled(f_cvc_pbc_minimum_image)) {
dist_v = main->center_of_mass() - ref1->center_of_mass();
} else {
dist_v = cvm::position_distance(ref1->center_of_mass(),
@ -244,14 +237,16 @@ void colvar::distance_z::calc_value()
}
} else {
if (b_no_PBC) {
if (!is_enabled(f_cvc_pbc_minimum_image)) {
dist_v = main->center_of_mass() -
(0.5 * (ref1->center_of_mass() + ref2->center_of_mass()));
axis = ref2->center_of_mass() - ref1->center_of_mass();
} else {
dist_v = cvm::position_distance(0.5 * (ref1->center_of_mass() +
ref2->center_of_mass()), main->center_of_mass());
axis = cvm::position_distance(ref1->center_of_mass(), ref2->center_of_mass());
ref2->center_of_mass()),
main->center_of_mass());
axis = cvm::position_distance(ref1->center_of_mass(),
ref2->center_of_mass());
}
axis_norm = axis.norm();
axis = axis.unit();
@ -268,16 +263,20 @@ void colvar::distance_z::calc_gradients()
if (fixed_axis) {
ref1->set_weighted_gradient(-1.0 * axis);
} else {
if (b_no_PBC) {
ref1->set_weighted_gradient( 1.0 / axis_norm * (main->center_of_mass() - ref2->center_of_mass() -
if (!is_enabled(f_cvc_pbc_minimum_image)) {
ref1->set_weighted_gradient( 1.0 / axis_norm *
(main->center_of_mass() - ref2->center_of_mass() -
x.real_value * axis ));
ref2->set_weighted_gradient( 1.0 / axis_norm * (ref1->center_of_mass() - main->center_of_mass() +
ref2->set_weighted_gradient( 1.0 / axis_norm *
(ref1->center_of_mass() - main->center_of_mass() +
x.real_value * axis ));
} else {
ref1->set_weighted_gradient( 1.0 / axis_norm * (
cvm::position_distance(ref2->center_of_mass(), main->center_of_mass()) - x.real_value * axis ));
cvm::position_distance(ref2->center_of_mass(),
main->center_of_mass()) - x.real_value * axis ));
ref2->set_weighted_gradient( 1.0 / axis_norm * (
cvm::position_distance(main->center_of_mass(), ref1->center_of_mass()) + x.real_value * axis ));
cvm::position_distance(main->center_of_mass(),
ref1->center_of_mass()) + x.real_value * axis ));
}
}
}
@ -390,17 +389,18 @@ colvar::distance_xy::distance_xy()
void colvar::distance_xy::calc_value()
{
if (b_no_PBC) {
if (!is_enabled(f_cvc_pbc_minimum_image)) {
dist_v = main->center_of_mass() - ref1->center_of_mass();
} else {
dist_v = cvm::position_distance(ref1->center_of_mass(),
main->center_of_mass());
}
if (!fixed_axis) {
if (b_no_PBC) {
if (!is_enabled(f_cvc_pbc_minimum_image)) {
v12 = ref2->center_of_mass() - ref1->center_of_mass();
} else {
v12 = cvm::position_distance(ref1->center_of_mass(), ref2->center_of_mass());
v12 = cvm::position_distance(ref1->center_of_mass(),
ref2->center_of_mass());
}
axis_norm = v12.norm();
axis = v12.unit();
@ -425,10 +425,11 @@ void colvar::distance_xy::calc_gradients()
ref1->set_weighted_gradient(-1.0 * x_inv * dist_v_ortho);
main->set_weighted_gradient( x_inv * dist_v_ortho);
} else {
if (b_no_PBC) {
if (!is_enabled(f_cvc_pbc_minimum_image)) {
v13 = main->center_of_mass() - ref1->center_of_mass();
} else {
v13 = cvm::position_distance(ref1->center_of_mass(), main->center_of_mass());
v13 = cvm::position_distance(ref1->center_of_mass(),
main->center_of_mass());
}
A = (dist_v * axis) / axis_norm;
@ -480,6 +481,7 @@ colvar::distance_dir::distance_dir(std::string const &conf)
{
function_type = "distance_dir";
enable(f_cvc_com_based);
enable(f_cvc_implicit_gradient);
x.type(colvarvalue::type_unit3vector);
}
@ -489,13 +491,14 @@ colvar::distance_dir::distance_dir()
{
function_type = "distance_dir";
enable(f_cvc_com_based);
enable(f_cvc_implicit_gradient);
x.type(colvarvalue::type_unit3vector);
}
void colvar::distance_dir::calc_value()
{
if (b_no_PBC) {
if (!is_enabled(f_cvc_pbc_minimum_image)) {
dist_v = group2->center_of_mass() - group1->center_of_mass();
} else {
dist_v = cvm::position_distance(group1->center_of_mass(),
@ -539,22 +542,26 @@ cvm::real colvar::distance_dir::dist2(colvarvalue const &x1,
colvarvalue colvar::distance_dir::dist2_lgrad(colvarvalue const &x1,
colvarvalue const &x2) const
{
return colvarvalue((x1.rvector_value - x2.rvector_value), colvarvalue::type_unit3vector);
return colvarvalue((x1.rvector_value - x2.rvector_value), colvarvalue::type_unit3vectorderiv);
}
colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const
{
return colvarvalue((x2.rvector_value - x1.rvector_value), colvarvalue::type_unit3vector);
return colvarvalue((x2.rvector_value - x1.rvector_value), colvarvalue::type_unit3vectorderiv);
}
colvar::distance_inv::distance_inv(std::string const &conf)
: distance(conf)
: cvc(conf)
{
function_type = "distance_inv";
group1 = parse_group(conf, "group1");
group2 = parse_group(conf, "group2");
get_keyval(conf, "exponent", exponent, 6);
if (exponent%2) {
cvm::error("Error: odd exponent provided, can only use even ones.\n");
@ -589,7 +596,7 @@ colvar::distance_inv::distance_inv()
void colvar::distance_inv::calc_value()
{
x.real_value = 0.0;
if (b_no_PBC) {
if (!is_enabled(f_cvc_pbc_minimum_image)) {
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
cvm::rvector const dv = ai2->pos - ai1->pos;
@ -655,14 +662,11 @@ colvar::distance_pairs::distance_pairs(std::string const &conf)
{
function_type = "distance_pairs";
if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) {
cvm::log("Computing distance using absolute positions (not minimal-image)");
}
group1 = parse_group(conf, "group1");
group2 = parse_group(conf, "group2");
x.type(colvarvalue::type_vector);
enable(f_cvc_implicit_gradient);
x.vector1d_value.resize(group1->size() * group2->size());
}
@ -670,6 +674,7 @@ colvar::distance_pairs::distance_pairs(std::string const &conf)
colvar::distance_pairs::distance_pairs()
{
function_type = "distance_pairs";
enable(f_cvc_implicit_gradient);
x.type(colvarvalue::type_vector);
}
@ -678,7 +683,7 @@ void colvar::distance_pairs::calc_value()
{
x.vector1d_value.resize(group1->size() * group2->size());
if (b_no_PBC) {
if (!is_enabled(f_cvc_pbc_minimum_image)) {
size_t i1, i2;
for (i1 = 0; i1 < group1->size(); i1++) {
for (i2 = 0; i2 < group2->size(); i2++) {
@ -693,7 +698,8 @@ void colvar::distance_pairs::calc_value()
size_t i1, i2;
for (i1 = 0; i1 < group1->size(); i1++) {
for (i2 = 0; i2 < group2->size(); i2++) {
cvm::rvector const dv = cvm::position_distance((*group1)[i1].pos, (*group2)[i2].pos);
cvm::rvector const dv = cvm::position_distance((*group1)[i1].pos,
(*group2)[i2].pos);
cvm::real const d = dv.norm();
x.vector1d_value[i1*group2->size() + i2] = d;
(*group1)[i1].grad = -1.0 * dv.unit();
@ -712,7 +718,7 @@ void colvar::distance_pairs::calc_gradients()
void colvar::distance_pairs::apply_force(colvarvalue const &force)
{
if (b_no_PBC) {
if (!is_enabled(f_cvc_pbc_minimum_image)) {
size_t i1, i2;
for (i1 = 0; i1 < group1->size(); i1++) {
for (i2 = 0; i2 < group2->size(); i2++) {
@ -725,7 +731,8 @@ void colvar::distance_pairs::apply_force(colvarvalue const &force)
size_t i1, i2;
for (i1 = 0; i1 < group1->size(); i1++) {
for (i2 = 0; i2 < group2->size(); i2++) {
cvm::rvector const dv = cvm::position_distance((*group1)[i1].pos, (*group2)[i2].pos);
cvm::rvector const dv = cvm::position_distance((*group1)[i1].pos,
(*group2)[i2].pos);
(*group1)[i1].apply_force(force[i1*group2->size() + i2] * (-1.0) * dv.unit());
(*group2)[i2].apply_force(force[i1*group2->size() + i2] * dv.unit());
}
@ -999,7 +1006,7 @@ colvar::rmsd::rmsd(std::string const &conf)
cvm::log("This is a standard minimum RMSD, derivatives of the optimal rotation "
"will not be computed as they cancel out in the gradients.");
atoms->b_fit_gradients = false;
atoms->disable(f_ag_fit_gradients);
// request the calculation of the derivatives of the rotation defined by the atom group
atoms->rot.request_group1_gradients(atoms->size());
@ -1191,7 +1198,7 @@ colvar::eigenvector::eigenvector(std::string const &conf)
atoms->b_rotate = true;
atoms->ref_pos = ref_pos;
atoms->center_ref_pos();
atoms->b_fit_gradients = false; // cancel out if group is fitted on itself
atoms->disable(f_ag_fit_gradients); // cancel out if group is fitted on itself
// and cvc is translationally invariant
// request the calculation of the derivatives of the rotation defined by the atom group
@ -1207,8 +1214,9 @@ colvar::eigenvector::eigenvector(std::string const &conf)
if (b_inline) {
cvm::log("Using vector components from input file.\n");
if (eigenvec.size() != atoms->size()) {
cvm::fatal_error("Error: vector components do not "
cvm::error("Error: vector components do not "
"match the number of requested atoms->\n");
return;
}
}
@ -1422,6 +1430,7 @@ colvar::cartesian::cartesian(std::string const &conf)
}
x.type(colvarvalue::type_vector);
enable(f_cvc_implicit_gradient);
x.vector1d_value.resize(atoms->size() * axes.size());
}

84
lib/colvars/colvarcomp_protein.cpp
View File

@ -20,15 +20,6 @@
// alpha component
//////////////////////////////////////////////////////////////////////
// FIXME: this will not make collect_gradients work
// because gradients in individual atom groups
// are those of the sub-cvcs (angle, hb), not those
// of this cvc (alpha)
// This is true of all cvcs with sub-cvcs, and those
// that do not calculate explicit gradients
// SO: we need a flag giving the availability of
// atomic gradients
colvar::alpha_angles::alpha_angles(std::string const &conf)
: cvc(conf)
{
@ -36,6 +27,7 @@ colvar::alpha_angles::alpha_angles(std::string const &conf)
cvm::log("Initializing alpha_angles object.\n");
function_type = "alpha_angles";
enable(f_cvc_implicit_gradient);
x.type(colvarvalue::type_scalar);
std::string segment_id;
@ -44,7 +36,7 @@ colvar::alpha_angles::alpha_angles(std::string const &conf)
std::vector<int> residues;
{
std::string residues_conf = "";
key_lookup(conf, "residueRange", residues_conf);
key_lookup(conf, "residueRange", &residues_conf);
if (residues_conf.size()) {
std::istringstream is(residues_conf);
int initial, final;
@ -57,12 +49,14 @@ colvar::alpha_angles::alpha_angles(std::string const &conf)
}
}
} else {
cvm::fatal_error("Error: no residues defined in \"residueRange\".\n");
cvm::error("Error: no residues defined in \"residueRange\".\n");
return;
}
}
if (residues.size() < 5) {
cvm::fatal_error("Error: not enough residues defined in \"residueRange\".\n");
cvm::error("Error: not enough residues defined in \"residueRange\".\n");
return;
}
std::string const &sid = segment_id;
@ -71,7 +65,8 @@ colvar::alpha_angles::alpha_angles(std::string const &conf)
get_keyval(conf, "hBondCoeff", hb_coeff, 0.5);
if ( (hb_coeff < 0.0) || (hb_coeff > 1.0) ) {
cvm::fatal_error("Error: hBondCoeff must be defined between 0 and 1.\n");
cvm::error("Error: hBondCoeff must be defined between 0 and 1.\n");
return;
}
@ -84,9 +79,9 @@ colvar::alpha_angles::alpha_angles(std::string const &conf)
theta.push_back(new colvar::angle(cvm::atom(r[i ], "CA", sid),
cvm::atom(r[i+1], "CA", sid),
cvm::atom(r[i+2], "CA", sid)));
atom_groups.push_back(theta.back()->atom_groups[0]);
atom_groups.push_back(theta.back()->atom_groups[1]);
atom_groups.push_back(theta.back()->atom_groups[2]);
register_atom_group(theta.back()->atom_groups[0]);
register_atom_group(theta.back()->atom_groups[1]);
register_atom_group(theta.back()->atom_groups[2]);
}
} else {
@ -106,7 +101,7 @@ colvar::alpha_angles::alpha_angles(std::string const &conf)
hb.push_back(new colvar::h_bond(cvm::atom(r[i ], "O", sid),
cvm::atom(r[i+4], "N", sid),
r0, en, ed));
atom_groups.push_back(hb.back()->atom_groups[0]);
register_atom_group(hb.back()->atom_groups[0]);
}
} else {
@ -123,6 +118,7 @@ colvar::alpha_angles::alpha_angles()
: cvc()
{
function_type = "alpha_angles";
enable(f_cvc_implicit_gradient);
x.type(colvarvalue::type_scalar);
}
@ -239,15 +235,6 @@ simple_scalar_dist_functions(alpha_angles)
// dihedral principal component
//////////////////////////////////////////////////////////////////////
// FIXME: this will not make collect_gradients work
// because gradients in individual atom groups
// are those of the sub-cvcs (dihedral), not those
// of this cvc
// This is true of all cvcs with sub-cvcs, and those
// that do not calculate explicit gradients
// SO: we need a flag giving the availability of
// atomic gradients
colvar::dihedPC::dihedPC(std::string const &conf)
: cvc(conf)
{
@ -255,6 +242,7 @@ colvar::dihedPC::dihedPC(std::string const &conf)
cvm::log("Initializing dihedral PC object.\n");
function_type = "dihedPC";
enable(f_cvc_implicit_gradient);
x.type(colvarvalue::type_scalar);
std::string segment_id;
@ -263,7 +251,7 @@ colvar::dihedPC::dihedPC(std::string const &conf)
std::vector<int> residues;
{
std::string residues_conf = "";
key_lookup(conf, "residueRange", residues_conf);
key_lookup(conf, "residueRange", &residues_conf);
if (residues_conf.size()) {
std::istringstream is(residues_conf);
int initial, final;
@ -276,12 +264,14 @@ colvar::dihedPC::dihedPC(std::string const &conf)
}
}
} else {
cvm::fatal_error("Error: no residues defined in \"residueRange\".\n");
cvm::error("Error: no residues defined in \"residueRange\".\n");
return;
}
}
if (residues.size() < 2) {
cvm::fatal_error("Error: dihedralPC requires at least two residues.\n");
cvm::error("Error: dihedralPC requires at least two residues.\n");
return;
}
std::string const &sid = segment_id;
@ -291,13 +281,16 @@ colvar::dihedPC::dihedPC(std::string const &conf)
int vecNumber;
if (get_keyval(conf, "vectorFile", vecFileName, vecFileName)) {
get_keyval(conf, "vectorNumber", vecNumber, 0);
if (vecNumber < 1)
cvm::fatal_error("A positive value of vectorNumber is required.");
if (vecNumber < 1) {
cvm::error("A positive value of vectorNumber is required.");
return;
}
std::ifstream vecFile;
vecFile.open(vecFileName.c_str());
if (!vecFile.good())
cvm::fatal_error("Error opening dihedral PCA vector file " + vecFileName + " for reading");
if (!vecFile.good()) {
cvm::error("Error opening dihedral PCA vector file " + vecFileName + " for reading");
}
// TODO: adapt to different formats by setting this flag
bool eigenvectors_as_columns = true;
@ -321,8 +314,9 @@ colvar::dihedPC::dihedPC(std::string const &conf)
for (int i = 1; i<vecNumber; i++)
vecFile.ignore(999999, '\n');
if (!vecFile.good())
cvm::fatal_error("Error reading dihedral PCA vector file " + vecFileName);
if (!vecFile.good()) {
cvm::error("Error reading dihedral PCA vector file " + vecFileName);
}
std::string line;
getline(vecFile, line);
@ -341,10 +335,11 @@ colvar::dihedPC::dihedPC(std::string const &conf)
}
if ( coeffs.size() != 4 * (residues.size() - 1)) {
cvm::fatal_error("Error: wrong number of coefficients: " +
cvm::error("Error: wrong number of coefficients: " +
cvm::to_str(coeffs.size()) + ". Expected " +
cvm::to_str(4 * (residues.size() - 1)) +
" (4 coeffs per residue, minus one residue).\n");
return;
}
for (size_t i = 0; i < residues.size()-1; i++) {
@ -353,19 +348,19 @@ colvar::dihedPC::dihedPC(std::string const &conf)
cvm::atom(r[i ], "CA", sid),
cvm::atom(r[i ], "C", sid),
cvm::atom(r[i+1], "N", sid)));
atom_groups.push_back(theta.back()->atom_groups[0]);
atom_groups.push_back(theta.back()->atom_groups[1]);
atom_groups.push_back(theta.back()->atom_groups[2]);
atom_groups.push_back(theta.back()->atom_groups[3]);
register_atom_group(theta.back()->atom_groups[0]);
register_atom_group(theta.back()->atom_groups[1]);
register_atom_group(theta.back()->atom_groups[2]);
register_atom_group(theta.back()->atom_groups[3]);
// Phi (next res)
theta.push_back(new colvar::dihedral(cvm::atom(r[i ], "C", sid),
cvm::atom(r[i+1], "N", sid),
cvm::atom(r[i+1], "CA", sid),
cvm::atom(r[i+1], "C", sid)));
atom_groups.push_back(theta.back()->atom_groups[0]);
atom_groups.push_back(theta.back()->atom_groups[1]);
atom_groups.push_back(theta.back()->atom_groups[2]);
atom_groups.push_back(theta.back()->atom_groups[3]);
register_atom_group(theta.back()->atom_groups[0]);
register_atom_group(theta.back()->atom_groups[1]);
register_atom_group(theta.back()->atom_groups[2]);
register_atom_group(theta.back()->atom_groups[3]);
}
if (cvm::debug())
@ -377,6 +372,7 @@ colvar::dihedPC::dihedPC()
: cvc()
{
function_type = "dihedPC";
enable(f_cvc_implicit_gradient);
x.type(colvarvalue::type_scalar);
}

14
lib/colvars/colvarcomp_rotations.cpp
View File

@ -22,6 +22,7 @@ colvar::orientation::orientation(std::string const &conf)
{
function_type = "orientation";
atoms = parse_group(conf, "atoms");
enable(f_cvc_implicit_gradient);
x.type(colvarvalue::type_quaternion);
ref_pos.reserve(atoms->size());
@ -29,8 +30,9 @@ colvar::orientation::orientation(std::string const &conf)
if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) {
cvm::log("Using reference positions from input file.\n");
if (ref_pos.size() != atoms->size()) {
cvm::fatal_error("Error: reference positions do not "
cvm::error("Error: reference positions do not "
"match the number of requested atoms.\n");
return;
}
}
@ -43,9 +45,11 @@ colvar::orientation::orientation(std::string const &conf)
if (get_keyval(conf, "refPositionsCol", file_col, std::string(""))) {
// use PDB flags if column is provided
bool found = get_keyval(conf, "refPositionsColValue", file_col_value, 0.0);
if (found && file_col_value==0.0)
cvm::fatal_error("Error: refPositionsColValue, "
if (found && file_col_value==0.0) {
cvm::error("Error: refPositionsColValue, "
"if provided, must be non-zero.\n");
return;
}
} else {
// if not, use atom indices
atoms->create_sorted_ids();
@ -56,8 +60,9 @@ colvar::orientation::orientation(std::string const &conf)
}
if (!ref_pos.size()) {
cvm::fatal_error("Error: must define a set of "
cvm::error("Error: must define a set of "
"reference coordinates.\n");
return;
}
@ -88,6 +93,7 @@ colvar::orientation::orientation()
: cvc()
{
function_type = "orientation";
enable(f_cvc_implicit_gradient);
x.type(colvarvalue::type_quaternion);
}

347
lib/colvars/colvardeps.cpp
View File

@ -10,24 +10,77 @@
#include "colvardeps.h"
colvardeps::colvardeps()
: time_step_factor (1) {}
colvardeps::~colvardeps() {
size_t i;
// Do not delete features if it's static
// for (i=0; i<features.size(); i++) {
// if (features[i] != NULL) delete features[i];
// }
remove_all_children();
// Protest if we are deleting an object while a parent object may still depend on it
// Another possible strategy is to have the child unlist itself from the parent's children
if (parents.size()) {
cvm::log("Warning: destroying " + description + " before its parents objects:");
cvm::log("Warning: destroying \"" + description + "\" before its parents objects:");
for (i=0; i<parents.size(); i++) {
cvm::log(parents[i]->description);
}
}
// Do not delete features if it's a static object
// may change in the future though
// for (i=0; i<features.size(); i++) {
// if (features[i] != NULL) delete features[i];
// }
remove_all_children();
}
void colvardeps::free_children_deps() {
// Dereference children requirements of all enabled features
// Useful when object is destroyed or set inactive
// CAUTION: when setting the parent object inactive, disable "active" first
// then call this, to avoid double-dereferencing the deps of "active"
// Cannot be in the base class destructor because it needs the derived class features()
size_t i,j,fid;
if (cvm::debug()) cvm::log("DEPS: freeing children deps for " + description);
cvm::increase_depth();
for (fid = 0; fid < feature_states.size(); fid++) {
if (is_enabled(fid)) {
for (i=0; i<features()[fid]->requires_children.size(); i++) {
int g = features()[fid]->requires_children[i];
for (j=0; j<children.size(); j++) {
if (cvm::debug()) cvm::log("DEPS: dereferencing children's "
+ children[j]->features()[g]->description);
children[j]->decr_ref_count(g);
}
}
}
}
cvm::decrease_depth();
}
// re-enable children features (and increase ref count accordingly)
// So free_children_deps() can be called whenever an object becomes inactive
void colvardeps::restore_children_deps() {
size_t i,j,fid;
cvm::increase_depth();
for (fid = 0; fid < feature_states.size(); fid++) {
if (is_enabled(fid)) {
for (i=0; i<features()[fid]->requires_children.size(); i++) {
int g = features()[fid]->requires_children[i];
for (j=0; j<children.size(); j++) {
if (cvm::debug()) cvm::log("DEPS: re-enabling children's "
+ children[j]->features()[g]->description);
children[j]->enable(g, false, false);
}
}
}
}
cvm::decrease_depth();
}
@ -37,15 +90,10 @@ void colvardeps::provide(int feature_id, bool truefalse) {
void colvardeps::set_enabled(int feature_id, bool truefalse) {
// if (!is_static(feature_id)) {
// cvm::error("Cannot set feature " + features()[feature_id]->description + " statically in " + description + ".\n");
// return;
// }
if (truefalse) {
// Resolve dependencies too
enable(feature_id);
} else {
feature_states[feature_id].enabled = false;
disable(feature_id);
}
}
@ -56,7 +104,7 @@ bool colvardeps::get_keyval_feature(colvarparse *cvp,
colvarparse::Parse_Mode const parse_mode)
{
if (!is_user(feature_id)) {
cvm::error("Cannot set feature " + features()[feature_id]->description + " from user input in " + description + ".\n");
cvm::error("Cannot set feature \"" + features()[feature_id]->description + "\" from user input in \"" + description + "\".\n");
return false;
}
bool value;
@ -83,21 +131,34 @@ int colvardeps::enable(int feature_id,
if (cvm::debug()) {
cvm::log("DEPS: " + description +
(dry_run ? " testing " : " requiring ") +
(dry_run ? " testing " : " enabling ") +
"\"" + f->description +"\"");
}
if (fs->enabled) {
// Do not try to solve deps if already enabled
if (!(dry_run || toplevel)) {
// This is a dependency
// Prevent disabling this feature as long
// as requirement is enabled
fs->ref_count++;
if (cvm::debug())
cvm::log("DEPS: bumping ref_count to " + cvm::to_str(fs->ref_count));
}
// Do not try to further resolve deps
return COLVARS_OK;
}
std::string feature_type_descr = is_static(feature_id) ? "Static" :
(is_dynamic(feature_id) ? "Dynamic" : "User-controlled");
if (!fs->available) {
if (!dry_run) {
if (toplevel) {
cvm::error("Error: Feature unavailable: \"" + f->description + "\" in " + description + ".");
cvm::error("Error: " + feature_type_descr + " feature unavailable: \""
+ f->description + "\" in " + description + ".");
} else {
cvm::log("Feature unavailable: \"" + f->description + "\" in " + description);
cvm::log(feature_type_descr + " feature unavailable: \""
+ f->description + "\" in " + description + ".");
}
}
return COLVARS_ERROR;
@ -105,21 +166,22 @@ int colvardeps::enable(int feature_id,
if (!toplevel && !is_dynamic(feature_id)) {
if (!dry_run) {
cvm::log("Non-dynamic feature : \"" + f->description
+ "\" in " + description + " may not be enabled as a dependency.\n");
cvm::log(feature_type_descr + " feature \"" + f->description
+ "\" may not be enabled as a dependency in " + description + ".\n");
}
return COLVARS_ERROR;
}
// 1) enforce exclusions
// reminder: exclusions must be mutual for this to work
for (i=0; i<f->requires_exclude.size(); i++) {
feature *g = features()[f->requires_exclude[i]];
if (cvm::debug())
cvm::log(f->description + " requires exclude " + g->description);
if (is_enabled(f->requires_exclude[i])) {
if (!dry_run) {
cvm::log("Features \"" + f->description + "\" is incompatible with \""
+ g->description + "\" in " + description);
cvm::log("Feature \"" + f->description + "\" is incompatible with \""
+ g->description + "\" in " + description + ".");
if (toplevel) {
cvm::error("Error: Failed dependency in " + description + ".");
}
@ -156,23 +218,27 @@ int colvardeps::enable(int feature_id,
res = enable(g, true, false); // see if available
if (res == COLVARS_OK) {
ok = true;
if (!dry_run) enable(g, false, false); // Require again, for real
if (!dry_run) {
enable(g, false, false); // Require again, for real
fs->alternate_refs.push_back(g); // We remember we enabled this
// so we can free it if this feature gets disabled
}
break;
}
}
if (!ok) {
if (!dry_run) {
cvm::log("No dependency satisfied among alternates:");
cvm::log("-----------------------------------------");
cvm::log("\"" + f->description + "\" in " + description
+ " requires one of the following features, none of which can be enabled:\n");
cvm::log("-----------------------------------------\n");
cvm::increase_depth();
for (j=0; j<f->requires_alt[i].size(); j++) {
int g = f->requires_alt[i][j];
cvm::log(cvm::to_str(j+1) + ". " + features()[g]->description);
cvm::increase_depth();
enable(g, false, false); // Just for printing error output
cvm::decrease_depth();
}
cvm::decrease_depth();
cvm::log("-----------------------------------------");
cvm::log("for \"" + f->description + "\" in " + description);
if (toplevel) {
cvm::error("Error: Failed dependency in " + description + ".");
}
@ -182,12 +248,13 @@ int colvardeps::enable(int feature_id,
}
// 4) solve deps in children
// if the object is inactive, we solve but do not enable: will be enabled
// when the object becomes active
cvm::increase_depth();
for (i=0; i<f->requires_children.size(); i++) {
int g = f->requires_children[i];
for (j=0; j<children.size(); j++) {
cvm::increase_depth();
res = children[j]->enable(g, dry_run, false);
cvm::decrease_depth();
res = children[j]->enable(g, dry_run || !is_enabled(), false);
if (res != COLVARS_OK) {
if (!dry_run) {
cvm::log("...required by \"" + f->description + "\" in " + description);
@ -198,25 +265,114 @@ int colvardeps::enable(int feature_id,
return res;
}
}
// If we've just touched the features of child objects, refresh them
if (!dry_run && f->requires_children.size() != 0) {
for (j=0; j<children.size(); j++) {
children[j]->refresh_deps();
}
}
}
cvm::decrease_depth();
// Actually enable feature only once everything checks out
if (!dry_run) fs->enabled = true;
if (!dry_run) {
fs->enabled = true;
// This should be the only reference
if (!toplevel) fs->ref_count = 1;
if (feature_id == 0) {
// Waking up this object, enable all deps in children
restore_children_deps();
}
do_feature_side_effects(feature_id);
if (cvm::debug())
cvm::log("DEPS: feature \"" + f->description + "\" in "
+ description + " enabled, ref_count = 1.");
}
return COLVARS_OK;
}
// disable() {
//
// // we need refs to parents to walk up the deps tree!
// // or refresh
// }
int colvardeps::disable(int feature_id) {
size_t i, j;
feature *f = features()[feature_id];
feature_state *fs = &feature_states[feature_id];
if (cvm::debug()) cvm::log("DEPS: disabling feature \""
+ f->description + "\" in " + description);
if (fs->enabled == false) {
return COLVARS_OK;
}
if (fs->ref_count > 1) {
cvm::error("Error: cannot disable feature \"" + f->description
+ "\" in " + description + " because of " + cvm::to_str(fs->ref_count-1)
+ " remaining references.\n" );
return COLVARS_ERROR;
}
// internal deps (self)
for (i=0; i<f->requires_self.size(); i++) {
if (cvm::debug()) cvm::log("DEPS: dereferencing self "
+ features()[f->requires_self[i]]->description);
decr_ref_count(f->requires_self[i]);
}
// alternates
for (i=0; i<fs->alternate_refs.size(); i++) {
if (cvm::debug()) cvm::log("DEPS: dereferencing alt "
+ features()[fs->alternate_refs[i]]->description);
decr_ref_count(fs->alternate_refs[i]);
}
// Forget these, now that they are dereferenced
fs->alternate_refs.clear();
// deps in children
// except if the object is inactive, then children dependencies
// have already been dereferenced by this function
// (or never referenced if feature was enabled while the object
// was inactive)
if (is_enabled()) {
cvm::increase_depth();
for (i=0; i<f->requires_children.size(); i++) {
int g = f->requires_children[i];
for (j=0; j<children.size(); j++) {
if (cvm::debug()) cvm::log("DEPS: dereferencing children's "
+ children[j]->features()[g]->description);
children[j]->decr_ref_count(g);
}
}
cvm::decrease_depth();
}
fs->enabled = false;
fs->ref_count = 0;
if (feature_id == 0) {
// Putting this object to sleep
free_children_deps();
}
return COLVARS_OK;
}
int colvardeps::decr_ref_count(int feature_id) {
int &rc = feature_states[feature_id].ref_count;
feature *f = features()[feature_id];
if (cvm::debug())
cvm::log("DEPS: decreasing reference count of \"" + f->description
+ "\" in " + description + ".\n");
if (rc <= 0) {
cvm::error("Error: cannot decrease reference count of feature \"" + f->description
+ "\" in " + description + ", which is " + cvm::to_str(rc) + ".\n");
return COLVARS_ERROR;
}
rc--;
if (rc == 0 && f->is_dynamic()) {
// we can auto-disable this feature
if (cvm::debug())
cvm::log("DEPS will now auto-disable dynamic feature \"" + f->description
+ "\" in " + description + ".\n");
disable(feature_id);
}
return COLVARS_OK;
}
void colvardeps::init_feature(int feature_id, const char *description, feature_type type) {
features()[feature_id]->description = description;
features()[feature_id]->type = type;
@ -235,6 +391,11 @@ void colvardeps::init_feature(int feature_id, const char *description, feature_t
features()[f]->requires_alt.back()[0] = g; \
features()[f]->requires_alt.back()[1] = h; \
features()[f]->requires_alt.back()[2] = i
#define f_req_alt4(f, g, h, i, j) features()[f]->requires_alt.push_back(std::vector<int>(4));\
features()[f]->requires_alt.back()[0] = g; \
features()[f]->requires_alt.back()[1] = h; \
features()[f]->requires_alt.back()[2] = i; \
features()[f]->requires_alt.back()[3] = j
void colvardeps::init_cvb_requires() {
int i;
@ -246,6 +407,9 @@ void colvardeps::init_cvb_requires() {
init_feature(f_cvb_active, "active", f_type_dynamic);
f_req_children(f_cvb_active, f_cv_active);
init_feature(f_cvb_awake, "awake", f_type_static);
f_req_self(f_cvb_awake, f_cvb_active);
init_feature(f_cvb_apply_force, "apply force", f_type_user);
f_req_children(f_cvb_apply_force, f_cv_gradient);
@ -278,9 +442,12 @@ void colvardeps::init_cv_requires() {
}
init_feature(f_cv_active, "active", f_type_dynamic);
f_req_children(f_cv_active, f_cvc_active);
// Colvars must be either a linear combination, or scalar (and polynomial) or scripted
f_req_alt3(f_cv_active, f_cv_scalar, f_cv_linear, f_cv_scripted);
// Do not require f_cvc_active in children, as some components may be disabled
// Colvars must be either a linear combination, or scalar (and polynomial) or scripted/custom
f_req_alt4(f_cv_active, f_cv_scalar, f_cv_linear, f_cv_scripted, f_cv_custom_function);
init_feature(f_cv_awake, "awake", f_type_static);
f_req_self(f_cv_awake, f_cv_active);
init_feature(f_cv_gradient, "gradient", f_type_dynamic);
f_req_children(f_cv_gradient, f_cvc_gradient);
@ -288,8 +455,10 @@ void colvardeps::init_cv_requires() {
init_feature(f_cv_collect_gradient, "collect gradient", f_type_dynamic);
f_req_self(f_cv_collect_gradient, f_cv_gradient);
f_req_self(f_cv_collect_gradient, f_cv_scalar);
// The following exlusion could be lifted by implementing the feature
f_req_exclude(f_cv_collect_gradient, f_cv_scripted);
init_feature(f_cv_fdiff_velocity, "fdiff_velocity", f_type_dynamic);
init_feature(f_cv_fdiff_velocity, "velocity from finite differences", f_type_dynamic);
// System force: either trivial (spring force); through extended Lagrangian, or calculated explicitly
init_feature(f_cv_total_force, "total force", f_type_dynamic);
@ -335,6 +504,9 @@ void colvardeps::init_cv_requires() {
init_feature(f_cv_subtract_applied_force, "subtract applied force from total force", f_type_user);
f_req_self(f_cv_subtract_applied_force, f_cv_total_force);
// There is no well-defined way to implement f_cv_subtract_applied_force
// in the case of extended-Lagrangian colvars
f_req_exclude(f_cv_subtract_applied_force, f_cv_extended_Lagrangian);
init_feature(f_cv_lower_boundary, "lower boundary", f_type_user);
f_req_self(f_cv_lower_boundary, f_cv_scalar);
@ -350,12 +522,21 @@ void colvardeps::init_cv_requires() {
init_feature(f_cv_corrfunc, "correlation function", f_type_user);
init_feature(f_cv_scripted, "scripted", f_type_static);
init_feature(f_cv_scripted, "scripted", f_type_user);
init_feature(f_cv_custom_function, "custom function", f_type_user);
f_req_exclude(f_cv_custom_function, f_cv_scripted);
init_feature(f_cv_periodic, "periodic", f_type_static);
f_req_self(f_cv_periodic, f_cv_homogeneous);
init_feature(f_cv_scalar, "scalar", f_type_static);
init_feature(f_cv_linear, "linear", f_type_static);
init_feature(f_cv_homogeneous, "homogeneous", f_type_static);
// because total forces are obtained from the previous time step,
// we cannot (currently) have colvar values and total forces for the same timestep
init_feature(f_cv_multiple_ts, "multiple timestep colvar");
f_req_exclude(f_cv_multiple_ts, f_cv_total_force_calc);
}
// Initialize feature_states for each instance
@ -365,23 +546,6 @@ void colvardeps::init_cv_requires() {
// Most features are available, so we set them so
// and list exceptions below
}
// // properties that may NOT be enabled as a dependency
// // This will be deprecated by feature types
// int unavailable_deps[] = {
// f_cv_lower_boundary,
// f_cv_upper_boundary,
// f_cv_extended_Lagrangian,
// f_cv_Langevin,
// f_cv_scripted,
// f_cv_periodic,
// f_cv_scalar,
// f_cv_linear,
// f_cv_homogeneous
// };
// for (i = 0; i < sizeof(unavailable_deps) / sizeof(unavailable_deps[0]); i++) {
// feature_states[unavailable_deps[i]].available = false;
// }
}
@ -401,20 +565,26 @@ void colvardeps::init_cvc_requires() {
init_feature(f_cvc_gradient, "gradient", f_type_dynamic);
init_feature(f_cvc_implicit_gradient, "implicit gradient", f_type_static);
f_req_children(f_cvc_implicit_gradient, f_ag_implicit_gradient);
init_feature(f_cvc_inv_gradient, "inverse gradient", f_type_dynamic);
f_req_self(f_cvc_inv_gradient, f_cvc_gradient);
init_feature(f_cvc_debug_gradient, "debug gradient", f_type_user);
f_req_self(f_cvc_debug_gradient, f_cvc_gradient);
f_req_exclude(f_cvc_debug_gradient, f_cvc_implicit_gradient);
init_feature(f_cvc_Jacobian, "Jacobian derivative", f_type_dynamic);
f_req_self(f_cvc_Jacobian, f_cvc_inv_gradient);
init_feature(f_cvc_com_based, "depends on group centers of mass", f_type_static);
// init_feature(f_cvc_pbc_minimum_image, "use minimum-image distances with PBCs", f_type_user);
// Compute total force on first site only to avoid unwanted
// coupling to other colvars (see e.g. Ciccotti et al., 2005)
init_feature(f_cvc_one_site_total_force, "compute total collective force only from one group center", f_type_user);
init_feature(f_cvc_one_site_total_force, "compute total force from one group", f_type_user);
f_req_self(f_cvc_one_site_total_force, f_cvc_com_based);
init_feature(f_cvc_scalable, "scalable calculation", f_type_static);
@ -438,11 +608,17 @@ void colvardeps::init_cvc_requires() {
feature_states.push_back(feature_state(avail, false));
}
// CVCs are enabled from the start - get disabled based on flags
feature_states[f_cvc_active].enabled = true;
// Features that are implemented by all cvcs by default
// Each cvc specifies what other features are available
feature_states[f_cvc_active].available = true;
feature_states[f_cvc_gradient].available = true;
// Use minimum-image distances by default
feature_states[f_cvc_pbc_minimum_image].enabled = true;
// Features that are implemented by default if their requirements are
feature_states[f_cvc_one_site_total_force].available = true;
@ -464,8 +640,10 @@ void colvardeps::init_ag_requires() {
init_feature(f_ag_center, "translational fit", f_type_static);
init_feature(f_ag_rotate, "rotational fit", f_type_static);
init_feature(f_ag_fitting_group, "reference positions group", f_type_static);
init_feature(f_ag_fit_gradient_group, "fit gradient for main group", f_type_static);
init_feature(f_ag_fit_gradient_ref, "fit gradient for reference group", f_type_static);
init_feature(f_ag_implicit_gradient, "implicit atom gradient", f_type_dynamic);
init_feature(f_ag_fit_gradients, "fit gradients", f_type_user);
f_req_exclude(f_ag_fit_gradients, f_ag_implicit_gradient);
init_feature(f_ag_atom_forces, "atomic forces", f_type_dynamic);
// parallel calculation implies that we have at least a scalable center of mass,
@ -493,29 +671,50 @@ void colvardeps::init_ag_requires() {
feature_states[f_ag_scalable_com].available = false;
// TODO make f_ag_scalable depend on f_ag_scalable_com (or something else)
feature_states[f_ag_scalable].available = true;
feature_states[f_ag_fit_gradients].available = true;
feature_states[f_ag_implicit_gradient].available = true;
}
void colvardeps::print_state() {
size_t i;
cvm::log("Enabled features of " + description);
cvm::log("Enabled features of \"" + description + "\" (with reference count)");
for (i = 0; i < feature_states.size(); i++) {
if (feature_states[i].enabled)
cvm::log("- " + features()[i]->description);
if (is_enabled(i))
cvm::log("- " + features()[i]->description + " ("
+ cvm::to_str(feature_states[i].ref_count) + ")");
}
cvm::increase_depth();
for (i=0; i<children.size(); i++) {
cvm::log("* child " + cvm::to_str(i+1));
cvm::increase_depth();
children[i]->print_state();
cvm::decrease_depth();
}
cvm::decrease_depth();
}
void colvardeps::add_child(colvardeps *child) {
children.push_back(child);
child->parents.push_back((colvardeps *)this);
// Solve dependencies of already enabled parent features
// in the new child
size_t i, fid;
cvm::increase_depth();
for (fid = 0; fid < feature_states.size(); fid++) {
if (is_enabled(fid)) {
for (i=0; i<features()[fid]->requires_children.size(); i++) {
int g = features()[fid]->requires_children[i];
if (cvm::debug()) cvm::log("DEPS: re-enabling children's "
+ child->features()[g]->description);
child->enable(g, false, false);
}
}
}
cvm::decrease_depth();
}

124
lib/colvars/colvardeps.h
View File

@ -23,10 +23,14 @@
/// 3. Static features are static properties of the object, determined
/// programatically at initialization time.
///
/// In all classes, feature 0 is active. When an object is inactivated
/// all its children dependencies are dereferenced (free_children_deps)
/// While the object is inactive, no dependency solving is done on children
/// it is done when the object is activated back (restore_children_deps)
class colvardeps {
public:
colvardeps() {}
colvardeps();
virtual ~colvardeps();
// Subclasses should initialize the following members:
@ -34,9 +38,10 @@ public:
std::string description; // reference to object name (cv, cvc etc.)
/// This contains the current state of each feature for each object
// since the feature class only contains static properties
struct feature_state {
feature_state(bool a, bool e)
: available(a), enabled(e) {}
: available(a), enabled(e), ref_count(0) {}
/// Feature may be enabled, subject to possible dependencies
bool available;
@ -44,9 +49,28 @@ public:
/// TODO consider implications for dependency solving: anyone who disables
/// it should trigger a refresh of parent objects
bool enabled; // see if this should be private depending on implementation
// bool enabledOnce; // this should trigger an update when object is evaluated
/// Number of features requiring this one as a dependency
/// When it falls to zero:
/// - a dynamic feature is disabled automatically
/// - other features may be disabled statically
int ref_count;
/// List of features that were enabled by this one
/// as part of an alternate requirement (for ref counting purposes)
/// This is necessary because we don't know which feature in the list
/// we enabled, otherwise
std::vector<int> alternate_refs;
};
protected:
/// Time step multiplier (for coarse-timestep biases & colvars)
/// Biases and colvars will only be calculated at those times
/// (f_cvb_awake and f_cv_awake); a
/// Biases use this to apply "impulse" biasing forces at the outer timestep
/// Unused by lower-level objects (cvcs and atom groups)
int time_step_factor;
private:
/// List of the states of all features
@ -61,10 +85,13 @@ private:
};
public:
/// \brief returns time_step_factor
inline int get_time_step_factor() const {return time_step_factor;}
/// Pair a numerical feature ID with a description and type
void init_feature(int feature_id, const char *description, feature_type type = f_type_not_set);
/// Describes a feature and its dependecies
/// Describes a feature and its dependencies
/// used in a static array within each subclass
class feature {
@ -120,30 +147,16 @@ public:
private:
// pointers to objects this object depends on
// list should be maintained by any code that modifies the object
// this could be secured by making lists of colvars / cvcs / atom groups private and modified through accessor functions
/// pointers to objects this object depends on
/// list should be maintained by any code that modifies the object
/// this could be secured by making lists of colvars / cvcs / atom groups private and modified through accessor functions
std::vector<colvardeps *> children;
// pointers to objects that depend on this object
// the size of this array is in effect a reference counter
/// pointers to objects that depend on this object
/// the size of this array is in effect a reference counter
std::vector<colvardeps *> parents;
public:
// disabling a feature f:
// if parents depend on f, tell them to refresh / check that they are ok?
// if children provide features to satisfy f ONLY, disable that
// When the state of this object has changed, recursively tell parents
// to enforce their dependencies
// void refresh_parents() {
//
// }
// std::vector<colvardeps *> parents; // Needed to trigger a refresh if capabilities of this object change
// End of members to be initialized by subclasses
// Checks whether given feature is enabled
// Defaults to querying f_*_active
inline bool is_enabled(int f = f_cv_active) const {
@ -161,9 +174,7 @@ public:
/// dependencies will be checked by enable()
void provide(int feature_id, bool truefalse = true);
/// Set the feature's enabled flag, without dependency check or resolution
/// To be used for static properties only
/// Checking for availability is up to the caller
/// Enable or disable, depending on flag value
void set_enabled(int feature_id, bool truefalse = true);
protected:
@ -178,31 +189,57 @@ protected:
public:
int enable(int f, bool dry_run = false, bool toplevel = true); // enable a feature and recursively solve its dependencies
// dry_run is set to true to recursively test if a feature is available, without enabling it
// int disable(int f);
/// enable a feature and recursively solve its dependencies
/// for proper reference counting, one should not add
/// spurious calls to enable()
/// dry_run is set to true to recursively test if a feature is available, without enabling it
int enable(int f, bool dry_run = false, bool toplevel = true);
/// Disable a feature, decrease the reference count of its dependencies
/// and recursively disable them as applicable
int disable(int f);
/// This function is called whenever feature states are changed outside
/// of the object's control, that is, by parents
/// Eventually it may also be used when properties of children change
virtual int refresh_deps() { return COLVARS_OK; }
/// disable all enabled features to free their dependencies
/// to be done when deleting the object
/// Cannot be in the base class destructor because it needs the derived class features()
void free_children_deps();
/// re-enable children features (to be used when object becomes active)
void restore_children_deps();
/// Decrement the reference count of a feature
/// disabling it if it's dynamic and count reaches zero
int decr_ref_count(int f);
/// Implements possible actions to be carried out
/// when a given feature is enabled
/// Base function does nothing, can be overloaded
virtual void do_feature_side_effects(int id) {}
// NOTE that all feature enums should start with f_*_active
enum features_biases {
/// \brief Bias is active
f_cvb_active,
f_cvb_apply_force, // will apply forces
f_cvb_get_total_force, // requires total forces
f_cvb_history_dependent, // depends on simulation history
f_cvb_scalar_variables, // requires scalar colvars
f_cvb_calc_pmf, // whether this bias will compute a PMF
/// \brief Bias is awake (active on its own accord) this timestep
f_cvb_awake,
/// \brief will apply forces
f_cvb_apply_force,
/// \brief requires total forces
f_cvb_get_total_force,
/// \brief depends on simulation history
f_cvb_history_dependent,
/// \brief requires scalar colvars
f_cvb_scalar_variables,
/// \brief whether this bias will compute a PMF
f_cvb_calc_pmf,
f_cvb_ntot
};
enum features_colvar {
/// \brief Calculate colvar
f_cv_active,
/// \brief Colvar is awake (active on its own accord) this timestep
f_cv_awake,
/// \brief Gradients are calculated and temporarily stored, so
/// that external forces can be applied
f_cv_gradient,
@ -254,12 +291,16 @@ public:
f_cv_corrfunc,
/// \brief Value and gradient computed by user script
f_cv_scripted,
/// \brief Value and gradient computed by user function through Lepton
f_cv_custom_function,
/// \brief Colvar is periodic
f_cv_periodic,
/// \brief is scalar
f_cv_scalar,
f_cv_linear,
f_cv_homogeneous,
/// \brief multiple timestep through time_step_factor
f_cv_multiple_ts,
/// \brief Number of colvar features
f_cv_ntot
};
@ -268,10 +309,13 @@ public:
f_cvc_active,
f_cvc_scalar,
f_cvc_gradient,
/// \brief CVC doesn't calculate and store explicit atom gradients
f_cvc_implicit_gradient,
f_cvc_inv_gradient,
/// \brief If enabled, calc_gradients() will call debug_gradients() for every group needed
f_cvc_debug_gradient,
f_cvc_Jacobian,
f_cvc_pbc_minimum_image,
f_cvc_one_site_total_force,
f_cvc_com_based,
f_cvc_scalable,
@ -287,9 +331,9 @@ public:
/// Perform a standard minimum msd fit for given atoms
/// ie. not using refpositionsgroup
// f_ag_min_msd_fit,
f_ag_fit_gradient_group,// TODO check that these are sometimes needed separately
// maybe for minimum RMSD?
f_ag_fit_gradient_ref,
/// \brief Does not have explicit atom gradients from parent CVC
f_ag_implicit_gradient,
f_ag_fit_gradients,
f_ag_atom_forces,
f_ag_scalable,
f_ag_scalable_com,

3
lib/colvars/colvargrid.cpp
View File

@ -144,7 +144,8 @@ void colvar_grid_gradient::write_1D_integral(std::ostream &os)
os << "# xi A(xi)\n";
if ( cv.size() != 1 ) {
cvm::fatal_error("Cannot write integral for multi-dimensional gradient grids.");
cvm::error("Cannot write integral for multi-dimensional gradient grids.");
return;
}
integral = 0.0;

6
lib/colvars/colvargrid.h
View File

@ -198,7 +198,6 @@ public:
/// Default constructor
colvar_grid() : has_data(false)
{
save_delimiters = false;
nd = nt = 0;
mult = 1;
this->setup();
@ -225,7 +224,6 @@ public:
widths(g.widths),
has_data(false)
{
save_delimiters = false;
}
/// \brief Constructor from explicit grid sizes \param nx_i Number
@ -237,7 +235,6 @@ public:
size_t mult_i = 1)
: has_data(false)
{
save_delimiters = false;
this->setup(nx_i, t, mult_i);
}
@ -248,7 +245,6 @@ public:
bool margin = false)
: has_data(false)
{
save_delimiters = false;
this->init_from_colvars(colvars, t, mult_i, margin);
}
@ -840,7 +836,7 @@ public:
// reallocate the array in case the grid params have just changed
if (new_params) {
init_from_boundaries();
// data.resize(0); // no longer needed: setup calls clear()
// data.clear(); // no longer needed: setup calls clear()
return this->setup(nx, T(), mult);
}

176
lib/colvars/colvarmodule.cpp
View File

@ -21,10 +21,14 @@
#include "colvarbias_meta.h"
#include "colvarbias_restraint.h"
#include "colvarscript.h"
#include "colvaratoms.h"
colvarmodule::colvarmodule(colvarproxy *proxy_in)
{
depth_s = 0;
cv_traj_os = NULL;
// pointer to the proxy object
if (proxy == NULL) {
proxy = proxy_in;
@ -33,12 +37,10 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
// TODO relax this error to handle multiple molecules in VMD
// once the module is not static anymore
cvm::error("Error: trying to allocate the collective "
"variable module twice.\n");
"variable module twice.\n", BUG_ERROR);
return;
}
depth_s = 0;
cvm::log(cvm::line_marker);
cvm::log("Initializing the collective variables module, version "+
cvm::to_str(COLVARS_VERSION)+".\n");
@ -222,9 +224,9 @@ int colvarmodule::parse_config(std::string &conf)
// update any necessary proxy data
proxy->setup();
if (cv_traj_os.is_open()) {
if (cv_traj_os != NULL) {
// configuration might have changed, better redo the labels
write_traj_label(cv_traj_os);
write_traj_label(*cv_traj_os);
}
return get_error();
@ -295,7 +297,7 @@ int colvarmodule::parse_colvars(std::string const &conf)
std::string colvar_conf = "";
size_t pos = 0;
while (parse->key_lookup(conf, "colvar", colvar_conf, pos)) {
while (parse->key_lookup(conf, "colvar", &colvar_conf, &pos)) {
if (colvar_conf.size()) {
cvm::log(cvm::line_marker);
@ -350,7 +352,7 @@ int colvarmodule::parse_biases_type(std::string const &conf,
{
std::string bias_conf = "";
size_t conf_saved_pos = 0;
while (parse->key_lookup(conf, keyword, bias_conf, conf_saved_pos)) {
while (parse->key_lookup(conf, keyword, &bias_conf, &conf_saved_pos)) {
if (bias_conf.size()) {
cvm::log(cvm::line_marker);
cvm::increase_depth();
@ -409,12 +411,6 @@ int colvarmodule::parse_biases(std::string const &conf)
size_t i;
for (i = 0; i < biases.size(); i++) {
biases[i]->enable(colvardeps::f_cvb_active);
if (cvm::debug())
biases[i]->print_state();
}
size_t n_hist_dep_biases = 0;
std::vector<std::string> hist_dep_biases_names;
for (i = 0; i < biases.size(); i++) {
@ -487,7 +483,8 @@ int colvarmodule::catch_input_errors(int result)
}
colvarbias * colvarmodule::bias_by_name(std::string const &name) {
colvarbias * colvarmodule::bias_by_name(std::string const &name)
{
colvarmodule *cv = cvm::main();
for (std::vector<colvarbias *>::iterator bi = cv->biases.begin();
bi != cv->biases.end();
@ -500,7 +497,8 @@ colvarbias * colvarmodule::bias_by_name(std::string const &name) {
}
colvar *colvarmodule::colvar_by_name(std::string const &name) {
colvar *colvarmodule::colvar_by_name(std::string const &name)
{
colvarmodule *cv = cvm::main();
for (std::vector<colvar *>::iterator cvi = cv->colvars.begin();
cvi != cv->colvars.end();
@ -513,6 +511,20 @@ colvar *colvarmodule::colvar_by_name(std::string const &name) {
}
cvm::atom_group *colvarmodule::atom_group_by_name(std::string const &name)
{
colvarmodule *cv = cvm::main();
for (std::vector<cvm::atom_group *>::iterator agi = cv->named_atom_groups.begin();
agi != cv->named_atom_groups.end();
agi++) {
if ((*agi)->name == name) {
return (*agi);
}
}
return NULL;
}
int colvarmodule::change_configuration(std::string const &bias_name,
std::string const &conf)
{
@ -521,7 +533,10 @@ int colvarmodule::change_configuration(std::string const &bias_name,
cvm::increase_depth();
colvarbias *b;
b = bias_by_name(bias_name);
if (b == NULL) { cvm::error("Error: bias not found: " + bias_name); }
if (b == NULL) {
cvm::error("Error: bias not found: " + bias_name);
return COLVARS_ERROR;
}
b->change_configuration(conf);
cvm::decrease_depth();
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
@ -534,7 +549,10 @@ std::string colvarmodule::read_colvar(std::string const &name)
colvar *c;
std::stringstream ss;
c = colvar_by_name(name);
if (c == NULL) { cvm::fatal_error("Error: colvar not found: " + name); }
if (c == NULL) {
cvm::error("Error: colvar not found: " + name);
return std::string();
}
ss << c->value();
cvm::decrease_depth();
return ss.str();
@ -547,7 +565,10 @@ cvm::real colvarmodule::energy_difference(std::string const &bias_name,
colvarbias *b;
cvm::real energy_diff = 0.;
b = bias_by_name(bias_name);
if (b == NULL) { cvm::fatal_error("Error: bias not found: " + bias_name); }
if (b == NULL) {
cvm::error("Error: bias not found: " + bias_name);
return 0.;
}
energy_diff = b->energy_difference(conf);
cvm::decrease_depth();
return energy_diff;
@ -666,27 +687,45 @@ int colvarmodule::calc_colvars()
cvm::log("Calculating collective variables.\n");
// calculate collective variables and their gradients
// First, we need to decide which biases are awake
// so they can activate colvars as needed
std::vector<colvarbias *>::iterator bi;
for (bi = biases.begin(); bi != biases.end(); bi++) {
int tsf = (*bi)->get_time_step_factor();
if (tsf > 0 && (step_absolute() % tsf == 0)) {
(*bi)->enable(colvardeps::f_cvb_awake);
} else {
(*bi)->disable(colvardeps::f_cvb_awake);
}
}
int error_code = COLVARS_OK;
std::vector<colvar *>::iterator cvi;
// Determine which colvars are active at this iteration
variables_active()->resize(0);
variables_active()->clear();
variables_active()->reserve(variables()->size());
for (cvi = variables()->begin(); cvi != variables()->end(); cvi++) {
// This is a dynamic feature - the next call should be to enable()
// or disable() when dynamic dependency resolution is fully implemented
(*cvi)->set_enabled(colvardeps::f_cv_active,
step_absolute() % (*cvi)->get_time_step_factor() == 0);
// Wake up or put to sleep variables
int tsf = (*cvi)->get_time_step_factor();
if (tsf > 0 && (step_absolute() % tsf == 0)) {
(*cvi)->enable(colvardeps::f_cv_awake);
} else {
(*cvi)->disable(colvardeps::f_cv_awake);
}
if ((*cvi)->is_enabled()) {
variables_active()->push_back(*cvi);
}
}
// if SMP support is available, split up the work
if (proxy->smp_enabled() == COLVARS_OK) {
// first, calculate how much work (currently, how many active CVCs) each colvar has
variables_active_smp()->resize(0);
variables_active_smp_items()->resize(0);
variables_active_smp()->clear();
variables_active_smp_items()->clear();
variables_active_smp()->reserve(variables_active()->size());
variables_active_smp_items()->reserve(variables_active()->size());
@ -748,7 +787,8 @@ int colvarmodule::calc_biases()
total_bias_energy = 0.0;
// update the list of active biases
biases_active()->resize(0);
// which may have changed based on f_cvb_awake in calc_colvars()
biases_active()->clear();
biases_active()->reserve(biases.size());
for (bi = biases.begin(); bi != biases.end(); bi++) {
if ((*bi)->is_enabled()) {
@ -828,8 +868,7 @@ int colvarmodule::update_colvar_forces()
"of colvars (if they have any).\n");
cvm::increase_depth();
for (cvi = variables()->begin(); cvi != variables()->end(); cvi++) {
// Here we call even inactive colvars, so they accumulate biasing forces
// as well as update their extended-system dynamics
// Inactive colvars will only reset their forces and return 0 energy
total_colvar_energy += (*cvi)->update_forces_energy();
if (cvm::get_error()) {
return COLVARS_ERROR;
@ -883,11 +922,13 @@ int colvarmodule::write_restart_files()
((cvm::step_absolute() % restart_out_freq) == 0) ) {
cvm::log("Writing the state file \""+
restart_out_name+"\".\n");
proxy->backup_file(restart_out_name.c_str());
restart_out_os.open(restart_out_name.c_str());
if (!restart_out_os.is_open() || !write_restart(restart_out_os))
cvm::error("Error: in writing restart file.\n");
restart_out_os.close();
proxy->backup_file(restart_out_name);
std::ostream *restart_out_os = proxy->output_stream(restart_out_name);
if (!restart_out_os) return cvm::get_error();
if (!write_restart(*restart_out_os)) {
return cvm::error("Error: in writing restart file.\n", FILE_ERROR);
}
proxy->close_output_stream(restart_out_name);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
@ -896,26 +937,26 @@ int colvarmodule::write_restart_files()
int colvarmodule::write_traj_files()
{
if (!cv_traj_os.is_open()) {
if (cv_traj_os == NULL) {
open_traj_file(cv_traj_name);
}
// write labels in the traj file every 1000 lines and at first timestep
if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || cvm::step_relative() == 0) {
write_traj_label(cv_traj_os);
write_traj_label(*cv_traj_os);
}
if ((cvm::step_absolute() % cv_traj_freq) == 0) {
write_traj(cv_traj_os);
write_traj(*cv_traj_os);
}
if (restart_out_freq && cv_traj_os.is_open()) {
if (restart_out_freq && (cv_traj_os != NULL)) {
// flush the trajectory file if we are at the restart frequency
if ( (cvm::step_relative() > 0) &&
((cvm::step_absolute() % restart_out_freq) == 0) ) {
cvm::log("Synchronizing (emptying the buffer of) trajectory file \""+
cv_traj_name+"\".\n");
cv_traj_os.flush();
proxy->flush_output_stream(cv_traj_os);
}
}
@ -1003,9 +1044,11 @@ int colvarmodule::reset()
index_groups.clear();
index_group_names.clear();
if (cv_traj_os.is_open()) {
proxy->reset();
if (cv_traj_os != NULL) {
// Do not close file here, as we might not be done with it yet.
cv_traj_os.flush();
proxy->flush_output_stream(cv_traj_os);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
@ -1264,9 +1307,9 @@ int colvarmodule::write_output_files()
}
cvm::decrease_depth();
if (cv_traj_os.is_open()) {
// do not close to avoid problems with multiple NAMD runs
cv_traj_os.flush();
if (cv_traj_os != NULL) {
// do not close, there may be another run command
proxy->flush_output_stream(cv_traj_os);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
@ -1380,9 +1423,10 @@ std::ostream & colvarmodule::write_restart(std::ostream &os)
return os;
}
int colvarmodule::open_traj_file(std::string const &file_name)
{
if (cv_traj_os.is_open()) {
if (cv_traj_os != NULL) {
return COLVARS_OK;
}
@ -1390,36 +1434,35 @@ int colvarmodule::open_traj_file(std::string const &file_name)
if (cv_traj_append) {
cvm::log("Appending to colvar trajectory file \""+file_name+
"\".\n");
cv_traj_os.open(file_name.c_str(), std::ios::app);
cv_traj_os = (cvm::proxy)->output_stream(file_name, std::ios::app);
} else {
cvm::log("Writing to colvar trajectory file \""+file_name+
"\".\n");
proxy->backup_file(file_name.c_str());
cv_traj_os.open(file_name.c_str());
cv_traj_os = (cvm::proxy)->output_stream(file_name);
}
if (!cv_traj_os.is_open()) {
if (cv_traj_os == NULL) {
cvm::error("Error: cannot write to file \""+file_name+"\".\n",
FILE_ERROR);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
return cvm::get_error();
}
int colvarmodule::close_traj_file()
{
if (cv_traj_os.is_open()) {
cv_traj_os.close();
if (cv_traj_os != NULL) {
proxy->close_output_stream(cv_traj_name);
cv_traj_os = NULL;
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
return cvm::get_error();
}
std::ostream & colvarmodule::write_traj_label(std::ostream &os)
{
if (!os.good()) {
cvm::error("Cannot write to trajectory file.");
return os;
}
os.setf(std::ios::scientific, std::ios::floatfield);
os << "# " << cvm::wrap_string("step", cvm::it_width-2)
@ -1437,13 +1480,16 @@ std::ostream & colvarmodule::write_traj_label(std::ostream &os)
(*bi)->write_traj_label(os);
}
os << "\n";
if (cvm::debug()) {
os.flush();
proxy->flush_output_stream(&os);
}
cvm::decrease_depth();
return os;
}
std::ostream & colvarmodule::write_traj(std::ostream &os)
{
os.setf(std::ios::scientific, std::ios::floatfield);
@ -1463,9 +1509,11 @@ std::ostream & colvarmodule::write_traj(std::ostream &os)
(*bi)->write_traj(os);
}
os << "\n";
if (cvm::debug()) {
os.flush();
proxy->flush_output_stream(&os);
}
cvm::decrease_depth();
return os;
}
@ -1540,25 +1588,19 @@ void colvarmodule::clear_error()
}
void cvm::error(std::string const &message, int code)
int colvarmodule::error(std::string const &message, int code)
{
set_error_bits(code);
proxy->error(message);
return get_error();
}
void cvm::fatal_error(std::string const &message)
int colvarmodule::fatal_error(std::string const &message)
{
// TODO once all non-fatal errors have been set to be handled by error(),
// set DELETE_COLVARS here for VMD to handle it
set_error_bits(FATAL_ERROR);
proxy->fatal_error(message);
}
void cvm::exit(std::string const &message)
{
proxy->exit(message);
return get_error();
}

62
lib/colvars/colvarmodule.h
View File

@ -10,9 +10,7 @@
#ifndef COLVARMODULE_H
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2017-03-09"
#endif
#include "colvars_version.h"
#ifndef COLVARS_DEBUG
#define COLVARS_DEBUG false
@ -54,11 +52,6 @@ You can browse the class hierarchy or the list of source files.
#include <vector>
#include <list>
#ifdef NAMD_VERSION
// use Lustre-friendly wrapper to POSIX write()
#include "fstream_namd.h"
#endif
class colvarparse;
class colvar;
class colvarbias;
@ -188,7 +181,13 @@ private:
/// Indexes of the items to calculate for each colvar
std::vector<int> colvars_smp_items;
/// Array of named atom groups
std::vector<atom_group *> named_atom_groups;
public:
/// Register a named atom group into named_atom_groups
inline void register_named_atom_group(atom_group * ag) {
named_atom_groups.push_back(ag);
}
/// Array of collective variables
std::vector<colvar *> *variables();
@ -319,12 +318,6 @@ public:
/// (Re)initialize the output trajectory and state file (does not write it yet)
int setup_output();
#ifdef NAMD_VERSION
typedef ofstream_namd ofstream;
#else
typedef std::ofstream ofstream;
#endif
/// Read the input restart file
std::istream & read_restart(std::istream &is);
/// Write the output restart file
@ -332,7 +325,7 @@ public:
/// Open a trajectory file if requested (and leave it open)
int open_traj_file(std::string const &file_name);
/// Close it
/// Close it (note: currently unused)
int close_traj_file();
/// Write in the trajectory file
std::ostream & write_traj(std::ostream &os);
@ -354,6 +347,9 @@ public:
/// Look up a colvar by name; returns NULL if not found
static colvar * colvar_by_name(std::string const &name);
/// Look up a named atom group by name; returns NULL if not found
static atom_group * atom_group_by_name(std::string const &name);
/// Load new configuration for the given bias -
/// currently works for harmonic (force constant and/or centers)
int change_configuration(std::string const &bias_name, std::string const &conf);
@ -452,10 +448,10 @@ public:
static void log(std::string const &message);
/// Print a message to the main log and exit with error code
static void fatal_error(std::string const &message);
static int fatal_error(std::string const &message);
/// Print a message to the main log and set global error code
static void error(std::string const &message, int code = COLVARS_ERROR);
static int error(std::string const &message, int code = COLVARS_ERROR);
/// Print a message to the main log and exit normally
static void exit(std::string const &message);
@ -473,7 +469,6 @@ public:
static rvector position_distance(atom_pos const &pos1,
atom_pos const &pos2);
/// \brief Get the square distance between two positions (with
/// periodic boundary conditions handled transparently)
///
@ -483,20 +478,6 @@ public:
static real position_dist2(atom_pos const &pos1,
atom_pos const &pos2);
/// \brief Get the closest periodic image to a reference position
/// \param pos The position to look for the closest periodic image
/// \param ref_pos (optional) The reference position
static void select_closest_image(atom_pos &pos,
atom_pos const &ref_pos);
/// \brief Perform select_closest_image() on a set of atomic positions
///
/// After that, distance vectors can then be calculated directly,
/// without using position_distance()
static void select_closest_images(std::vector<atom_pos> &pos,
atom_pos const &ref_pos);
/// \brief Names of groups from a Gromacs .ndx file to be read at startup
std::list<std::string> index_group_names;
@ -556,14 +537,11 @@ protected:
std::string cv_traj_name;
/// Collective variables output trajectory file
colvarmodule::ofstream cv_traj_os;
std::ostream *cv_traj_os;
/// Appending to the existing trajectory file?
bool cv_traj_append;
/// Output restart file
colvarmodule::ofstream restart_out_os;
private:
/// Counter for the current depth in the object hierarchy (useg e.g. in output)
@ -704,18 +682,6 @@ inline void cvm::request_total_force()
proxy->request_total_force(true);
}
inline void cvm::select_closest_image(atom_pos &pos,
atom_pos const &ref_pos)
{
proxy->select_closest_image(pos, ref_pos);
}
inline void cvm::select_closest_images(std::vector<atom_pos> &pos,
atom_pos const &ref_pos)
{
proxy->select_closest_images(pos, ref_pos);
}
inline cvm::rvector cvm::position_distance(atom_pos const &pos1,
atom_pos const &pos2)
{

49
lib/colvars/colvarparse.cpp
View File

@ -17,10 +17,7 @@
// space & tab
std::string const colvarparse::white_space = " \t";
std::string colvarparse::dummy_string = "";
size_t colvarparse::dummy_pos = 0;
char const * const colvarparse::white_space = " \t";
// definition of single-value keyword parsers
@ -37,7 +34,7 @@ template<typename TYPE> bool colvarparse::_get_keyval_scalar_(std::string const
do {
std::string data_this = "";
b_found = key_lookup(conf, key, data_this, save_pos);
b_found = key_lookup(conf, key, &data_this, &save_pos);
if (b_found) {
if (!b_found_any)
b_found_any = true;
@ -92,7 +89,7 @@ bool colvarparse::_get_keyval_scalar_string_(std::string const &conf,
do {
std::string data_this = "";
b_found = key_lookup(conf, key, data_this, save_pos);
b_found = key_lookup(conf, key, &data_this, &save_pos);
if (b_found) {
if (!b_found_any)
b_found_any = true;
@ -156,7 +153,7 @@ template<typename TYPE> bool colvarparse::_get_keyval_vector_(std::string const
do {
std::string data_this = "";
b_found = key_lookup(conf, key, data_this, save_pos);
b_found = key_lookup(conf, key, &data_this, &save_pos);
if (b_found) {
if (!b_found_any)
b_found_any = true;
@ -313,7 +310,7 @@ bool colvarparse::get_keyval(std::string const &conf,
do {
std::string data_this = "";
b_found = key_lookup(conf, key, data_this, save_pos);
b_found = key_lookup(conf, key, &data_this, &save_pos);
if (b_found) {
if (!b_found_any)
b_found_any = true;
@ -552,8 +549,8 @@ std::istream & colvarparse::getline_nocomments(std::istream &is,
bool colvarparse::key_lookup(std::string const &conf,
char const *key_in,
std::string &data,
size_t &save_pos)
std::string *data,
size_t *save_pos)
{
if (cvm::debug()) {
cvm::log("Looking for the keyword \""+std::string(key_in)+"\" and its value.\n");
@ -570,14 +567,12 @@ bool colvarparse::key_lookup(std::string const &conf,
std::string const conf_lower(to_lower_cppstr(conf));
// by default, there is no value, unless we found one
data = "";
// when the function is invoked without save_pos, ensure that we
// start from zero
colvarparse::dummy_pos = 0;
if (data != NULL) {
data->clear();
}
// start from the first occurrence of key
size_t pos = conf_lower.find(key, save_pos);
size_t pos = conf_lower.find(key, (save_pos != NULL) ? *save_pos : 0);
// iterate over all instances of the substring until it finds it as isolated keyword
while (true) {
@ -593,7 +588,7 @@ bool colvarparse::key_lookup(std::string const &conf,
bool b_isolated_left = true, b_isolated_right = true;
if (pos > 0) {
if ( std::string("\n"+white_space+
if ( std::string("\n"+std::string(white_space)+
"}").find(conf[pos-1]) ==
std::string::npos ) {
// none of the valid delimiting characters is on the left of key
@ -602,7 +597,7 @@ bool colvarparse::key_lookup(std::string const &conf,
}
if (pos < conf.size()-key.size()-1) {
if ( std::string("\n"+white_space+
if ( std::string("\n"+std::string(white_space)+
"{").find(conf[pos+key.size()]) ==
std::string::npos ) {
// none of the valid delimiting characters is on the right of key
@ -625,9 +620,11 @@ bool colvarparse::key_lookup(std::string const &conf,
}
}
if (save_pos != NULL) {
// save the pointer for a future call (when iterating over multiple
// valid instances of the same keyword)
save_pos = pos + key.size();
*save_pos = pos + key.size();
}
// get the remainder of the line
size_t pl = conf.rfind("\n", pos);
@ -716,19 +713,21 @@ bool colvarparse::key_lookup(std::string const &conf,
data_end) + 1;
}
data.append(line, data_begin, (data_end-data_begin));
if (data != NULL) {
data->append(line, data_begin, (data_end-data_begin));
if (cvm::debug()) {
cvm::log("Keyword value = \""+data+"\".\n");
cvm::log("Keyword value = \""+*data+"\".\n");
}
if (data.size() && save_delimiters) {
data_begin_pos.push_back(conf.find(data, pos+key.size()));
data_end_pos.push_back(data_begin_pos.back()+data.size());
if (data->size()) {
data_begin_pos.push_back(conf.find(*data, pos+key.size()));
data_end_pos.push_back(data_begin_pos.back()+data->size());
}
}
}
save_pos = line_end;
if (save_pos != NULL) *save_pos = line_end;
return true;
}

25
lib/colvars/colvarparse.h
View File

@ -24,7 +24,7 @@
/// need to parse input inherit from this
class colvarparse {
protected:
private:
/// \brief List of legal keywords for this object: this is updated
/// by each call to colvarparse::get_keyval() or
@ -41,14 +41,6 @@ protected:
/// values before the keyword check is performed
std::list<size_t> data_end_pos;
/// \brief Whether or not to accumulate data_begin_pos and
/// data_end_pos in key_lookup(); it may be useful to disable
/// this after the constructor is called, because other files may be
/// read (e.g. restart) that would mess up the registry; in any
/// case, nothing serious happens until check_keywords() is invoked
/// (which should happen only right after construction)
bool save_delimiters;
/// \brief Add a new valid keyword to the list
void add_keyword(char const *key);
@ -62,14 +54,12 @@ public:
inline colvarparse()
: save_delimiters(true)
{
init();
}
/// Constructor that stores the object's config string
inline colvarparse(const std::string& conf)
: save_delimiters(true)
{
init(conf);
}
@ -115,8 +105,6 @@ public:
/// \brief Use this after parsing a config string (note that check_keywords() calls it already)
void clear_keyword_registry();
public:
/// \fn get_keyval bool const get_keyval (std::string const &conf,
/// char const *key, _type_ &value, _type_ const &def_value,
/// Parse_Mode const parse_mode) \brief Helper function to parse
@ -282,7 +270,7 @@ public:
/// Accepted white space delimiters, used in key_lookup()
static std::string const white_space;
static const char * const white_space;
/// \brief Low-level function for parsing configuration strings;
/// automatically adds the requested keyword to the list of valid
@ -293,13 +281,8 @@ public:
/// within "conf", useful when doing multiple calls
bool key_lookup(std::string const &conf,
char const *key,
std::string &data = dummy_string,
size_t &save_pos = dummy_pos);
/// Used as a default argument by key_lookup
static std::string dummy_string;
/// Used as a default argument by key_lookup
static size_t dummy_pos;
std::string *data = NULL,
size_t *save_pos = NULL);
/// \brief Works as std::getline() but also removes everything
/// between a comment character and the following newline

492
lib/colvars/colvarproxy.cpp Normal file
View File

@ -0,0 +1,492 @@
// -*- c++ -*-
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/colvars/colvars
// Please update all Colvars source files before making any changes.
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#include <sstream>
#include <string.h>
#include "colvarmodule.h"
#include "colvarproxy.h"
#include "colvarscript.h"
#include "colvaratoms.h"
colvarproxy_system::colvarproxy_system() {}
colvarproxy_system::~colvarproxy_system() {}
void colvarproxy_system::add_energy(cvm::real energy) {}
void colvarproxy_system::request_total_force(bool yesno)
{
if (yesno == true)
cvm::error("Error: total forces are currently not implemented.\n",
COLVARS_NOT_IMPLEMENTED);
}
bool colvarproxy_system::total_forces_enabled() const
{
return false;
}
cvm::real colvarproxy_system::position_dist2(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2)
{
return (position_distance(pos1, pos2)).norm2();
}
colvarproxy_atoms::colvarproxy_atoms() {}
colvarproxy_atoms::~colvarproxy_atoms()
{
reset();
}
int colvarproxy_atoms::reset()
{
atoms_ids.clear();
atoms_ncopies.clear();
atoms_masses.clear();
atoms_charges.clear();
atoms_positions.clear();
atoms_total_forces.clear();
atoms_new_colvar_forces.clear();
return COLVARS_OK;
}
int colvarproxy_atoms::add_atom_slot(int atom_id)
{
atoms_ids.push_back(atom_id);
atoms_ncopies.push_back(1);
atoms_masses.push_back(1.0);
atoms_charges.push_back(0.0);
atoms_positions.push_back(cvm::rvector(0.0, 0.0, 0.0));
atoms_total_forces.push_back(cvm::rvector(0.0, 0.0, 0.0));
atoms_new_colvar_forces.push_back(cvm::rvector(0.0, 0.0, 0.0));
return (atoms_ids.size() - 1);
}
int colvarproxy_atoms::init_atom(cvm::residue_id const &residue,
std::string const &atom_name,
std::string const &segment_id)
{
cvm::error("Error: initializing an atom by name and residue number is currently not supported.\n",
COLVARS_NOT_IMPLEMENTED);
return COLVARS_NOT_IMPLEMENTED;
}
int colvarproxy_atoms::check_atom_id(cvm::residue_id const &residue,
std::string const &atom_name,
std::string const &segment_id)
{
colvarproxy_atoms::init_atom(residue, atom_name, segment_id);
return COLVARS_NOT_IMPLEMENTED;
}
void colvarproxy_atoms::clear_atom(int index)
{
if (((size_t) index) >= atoms_ids.size()) {
cvm::error("Error: trying to disable an atom that was not previously requested.\n",
INPUT_ERROR);
}
if (atoms_ncopies[index] > 0) {
atoms_ncopies[index] -= 1;
}
}
int colvarproxy_atoms::load_atoms(char const *filename,
cvm::atom_group &atoms,
std::string const &pdb_field,
double const)
{
return cvm::error("Error: loading atom identifiers from a file "
"is currently not implemented.\n",
COLVARS_NOT_IMPLEMENTED);
}
int colvarproxy_atoms::load_coords(char const *filename,
std::vector<cvm::atom_pos> &pos,
const std::vector<int> &indices,
std::string const &pdb_field,
double const)
{
return cvm::error("Error: loading atomic coordinates from a file "
"is currently not implemented.\n",
COLVARS_NOT_IMPLEMENTED);
}
colvarproxy_atom_groups::colvarproxy_atom_groups() {}
colvarproxy_atom_groups::~colvarproxy_atom_groups()
{
reset();
}
int colvarproxy_atom_groups::reset()
{
atom_groups_ids.clear();
atom_groups_ncopies.clear();
atom_groups_masses.clear();
atom_groups_charges.clear();
atom_groups_coms.clear();
atom_groups_total_forces.clear();
atom_groups_new_colvar_forces.clear();
return COLVARS_OK;
}
int colvarproxy_atom_groups::add_atom_group_slot(int atom_group_id)
{
atom_groups_ids.push_back(atom_group_id);
atom_groups_ncopies.push_back(1);
atom_groups_masses.push_back(1.0);
atom_groups_charges.push_back(0.0);
atom_groups_coms.push_back(cvm::rvector(0.0, 0.0, 0.0));
atom_groups_total_forces.push_back(cvm::rvector(0.0, 0.0, 0.0));
atom_groups_new_colvar_forces.push_back(cvm::rvector(0.0, 0.0, 0.0));
return (atom_groups_ids.size() - 1);
}
int colvarproxy_atom_groups::scalable_group_coms()
{
return COLVARS_NOT_IMPLEMENTED;
}
int colvarproxy_atom_groups::init_atom_group(std::vector<int> const &atoms_ids)
{
cvm::error("Error: initializing a group outside of the Colvars module "
"is currently not supported.\n",
COLVARS_NOT_IMPLEMENTED);
return COLVARS_NOT_IMPLEMENTED;
}
void colvarproxy_atom_groups::clear_atom_group(int index)
{
if (((size_t) index) >= atom_groups_ids.size()) {
cvm::error("Error: trying to disable an atom group "
"that was not previously requested.\n",
INPUT_ERROR);
}
if (atom_groups_ncopies[index] > 0) {
atom_groups_ncopies[index] -= 1;
}
}
colvarproxy_smp::colvarproxy_smp()
{
b_smp_active = true;
}
colvarproxy_smp::~colvarproxy_smp() {}
int colvarproxy_smp::smp_enabled()
{
return COLVARS_NOT_IMPLEMENTED;
}
int colvarproxy_smp::smp_colvars_loop()
{
return COLVARS_NOT_IMPLEMENTED;
}
int colvarproxy_smp::smp_biases_loop()
{
return COLVARS_NOT_IMPLEMENTED;
}
int colvarproxy_smp::smp_biases_script_loop()
{
return COLVARS_NOT_IMPLEMENTED;
}
int colvarproxy_smp::smp_thread_id()
{
return COLVARS_NOT_IMPLEMENTED;
}
int colvarproxy_smp::smp_num_threads()
{
return COLVARS_NOT_IMPLEMENTED;
}
int colvarproxy_smp::smp_lock()
{
return COLVARS_OK;
}
int colvarproxy_smp::smp_trylock()
{
return COLVARS_OK;
}
int colvarproxy_smp::smp_unlock()
{
return COLVARS_OK;
}
colvarproxy_replicas::colvarproxy_replicas() {}
colvarproxy_replicas::~colvarproxy_replicas() {}
bool colvarproxy_replicas::replica_enabled()
{
return false;
}
int colvarproxy_replicas::replica_index()
{
return 0;
}
int colvarproxy_replicas::replica_num()
{
return 1;
}
void colvarproxy_replicas::replica_comm_barrier() {}
int colvarproxy_replicas::replica_comm_recv(char* msg_data,
int buf_len,
int src_rep)
{
return COLVARS_NOT_IMPLEMENTED;
}
int colvarproxy_replicas::replica_comm_send(char* msg_data,
int msg_len,
int dest_rep)
{
return COLVARS_NOT_IMPLEMENTED;
}
colvarproxy_script::colvarproxy_script()
{
script = NULL;
}
colvarproxy_script::~colvarproxy_script() {}
char *colvarproxy_script::script_obj_to_str(unsigned char *obj)
{
return reinterpret_cast<char *>(obj);
}
int colvarproxy_script::run_force_callback()
{
return COLVARS_NOT_IMPLEMENTED;
}
int colvarproxy_script::run_colvar_callback(
std::string const &name,
std::vector<const colvarvalue *> const &cvcs,
colvarvalue &value)
{
return COLVARS_NOT_IMPLEMENTED;
}
int colvarproxy_script::run_colvar_gradient_callback(
std::string const &name,
std::vector<const colvarvalue *> const &cvcs,
std::vector<cvm::matrix2d<cvm::real> > &gradient)
{
return COLVARS_NOT_IMPLEMENTED;
}
colvarproxy_io::colvarproxy_io() {}
colvarproxy_io::~colvarproxy_io() {}
int colvarproxy_io::get_frame(long int&)
{
return COLVARS_NOT_IMPLEMENTED;
}
int colvarproxy_io::set_frame(long int)
{
return COLVARS_NOT_IMPLEMENTED;
}
std::ostream * colvarproxy_io::output_stream(std::string const &output_name,
std::ios_base::openmode mode)
{
if (cvm::debug()) {
cvm::log("Using colvarproxy::output_stream()\n");
}
std::list<std::ostream *>::iterator osi = output_files.begin();
std::list<std::string>::iterator osni = output_stream_names.begin();
for ( ; osi != output_files.end(); osi++, osni++) {
if (*osni == output_name) {
return *osi;
}
}
if (!(mode & (std::ios_base::app | std::ios_base::ate))) {
backup_file(output_name);
}
std::ofstream *os = new std::ofstream(output_name.c_str(), mode);
if (!os->is_open()) {
cvm::error("Error: cannot write to file/channel \""+output_name+"\".\n",
FILE_ERROR);
return NULL;
}
output_stream_names.push_back(output_name);
output_files.push_back(os);
return os;
}
int colvarproxy_io::flush_output_stream(std::ostream *os)
{
std::list<std::ostream *>::iterator osi = output_files.begin();
std::list<std::string>::iterator osni = output_stream_names.begin();
for ( ; osi != output_files.end(); osi++, osni++) {
if (*osi == os) {
((std::ofstream *) (*osi))->flush();
return COLVARS_OK;
}
}
return cvm::error("Error: trying to flush an output file/channel "
"that wasn't open.\n", BUG_ERROR);
}
int colvarproxy_io::close_output_stream(std::string const &output_name)
{
std::list<std::ostream *>::iterator osi = output_files.begin();
std::list<std::string>::iterator osni = output_stream_names.begin();
for ( ; osi != output_files.end(); osi++, osni++) {
if (*osni == output_name) {
((std::ofstream *) (*osi))->close();
output_files.erase(osi);
output_stream_names.erase(osni);
return COLVARS_OK;
}
}
return cvm::error("Error: trying to close an output file/channel "
"that wasn't open.\n", BUG_ERROR);
}
int colvarproxy_io::backup_file(char const *filename)
{
return COLVARS_NOT_IMPLEMENTED;
}
colvarproxy::colvarproxy()
{
colvars = NULL;
b_simulation_running = true;
}
colvarproxy::~colvarproxy() {}
int colvarproxy::reset()
{
int error_code = COLVARS_OK;
error_code |= colvarproxy_atoms::reset();
error_code |= colvarproxy_atom_groups::reset();
return error_code;
}
int colvarproxy::setup()
{
return COLVARS_OK;
}
int colvarproxy::update_input()
{
return COLVARS_OK;
}
int colvarproxy::update_output()
{
return COLVARS_OK;
}
size_t colvarproxy::restart_frequency()
{
return 0;
}

868
lib/colvars/colvarproxy.h
View File

@ -16,55 +16,36 @@
#include "colvarmodule.h"
#include "colvarvalue.h"
/// \file colvarproxy.h
/// \brief Colvars proxy classes
///
/// This file declares the class for the object responsible for interfacing
/// Colvars with other codes (MD engines, VMD, Python). The \link colvarproxy
/// \endlink class is a derivative of multiple classes, each devoted to a
/// specific task (e.g. \link colvarproxy_atoms \endlink to access data for
/// individual atoms).
///
/// To interface to a new MD engine, the simplest solution is to derive a new
/// class from \link colvarproxy \endlink. Currently implemented are: \link
/// colvarproxy_lammps, \endlink, \link colvarproxy_namd, \endlink, \link
/// colvarproxy_vmd, \endlink.
// forward declarations
class colvarscript;
/// \brief Interface between the collective variables module and
/// the simulation or analysis program (NAMD, VMD, LAMMPS...).
/// This is the base class: each interfaced program is supported by a derived class.
/// Only pure virtual functions ("= 0") must be reimplemented to ensure baseline functionality.
class colvarproxy {
/// Methods for accessing the simulation system (PBCs, integrator, etc)
class colvarproxy_system {
public:
/// Pointer to the main object
colvarmodule *colvars;
/// Constructor
colvarproxy()
{
colvars = NULL;
b_simulation_running = true;
b_smp_active = true;
script = NULL;
}
colvarproxy_system();
/// Destructor
virtual ~colvarproxy()
{}
/// (Re)initialize required member data after construction
virtual int setup()
{
return COLVARS_OK;
}
/// \brief Update data required by the colvars module (e.g. cache atom positions)
///
/// TODO Break up colvarproxy_namd and colvarproxy_lammps function into these
virtual int update_input()
{
return COLVARS_OK;
}
/// \brief Update data based from the results of a module update (e.g. send forces)
virtual int update_output()
{
return COLVARS_OK;
}
// **************** SIMULATION PARAMETERS ****************
virtual ~colvarproxy_system();
/// \brief Value of the unit for atomic coordinates with respect to
/// angstroms (used by some variables for hard-coded default values)
@ -73,7 +54,7 @@ public:
/// \brief Boltzmann constant
virtual cvm::real boltzmann() = 0;
/// \brief Temperature of the simulation (K)
/// \brief Target temperature of the simulation (K units)
virtual cvm::real temperature() = 0;
/// \brief Time step of the simulation (fs)
@ -82,263 +63,9 @@ public:
/// \brief Pseudo-random number with Gaussian distribution
virtual cvm::real rand_gaussian(void) = 0;
/// \brief Get the current frame number
// Returns error code
virtual int get_frame(long int&) { return COLVARS_NOT_IMPLEMENTED; }
/// \brief Set the current frame number (as well as colvarmodule::it)
// Returns error code
virtual int set_frame(long int) { return COLVARS_NOT_IMPLEMENTED; }
/// \brief Prefix to be used for input files (restarts, not
/// configuration)
std::string input_prefix_str, output_prefix_str, restart_output_prefix_str;
inline std::string & input_prefix()
{
return input_prefix_str;
}
/// \brief Prefix to be used for output restart files
inline std::string restart_output_prefix()
{
return restart_output_prefix_str;
}
/// \brief Prefix to be used for output files (final system
/// configuration)
inline std::string output_prefix()
{
return output_prefix_str;
}
/// \brief Restarts will be written each time this number of steps has passed
virtual size_t restart_frequency()
{
return 0;
}
protected:
/// Whether a simulation is running (and try to prevent irrecovarable errors)
bool b_simulation_running;
public:
/// Whether a simulation is running (and try to prevent irrecovarable errors)
virtual bool simulation_running() const
{
return b_simulation_running;
}
protected:
/// \brief Currently opened output files: by default, these are ofstream objects.
/// Allows redefinition to implement different output mechanisms
std::list<std::ostream *> output_files;
/// \brief Identifiers for output_stream objects: by default, these are the names of the files
std::list<std::string> output_stream_names;
public:
// ***************** SHARED-MEMORY PARALLELIZATION *****************
/// Whether threaded parallelization is available (TODO: make this a cvm::deps feature)
virtual int smp_enabled()
{
return COLVARS_NOT_IMPLEMENTED;
}
/// Whether threaded parallelization should be used (TODO: make this a cvm::deps feature)
bool b_smp_active;
/// Distribute calculation of colvars (and their components) across threads
virtual int smp_colvars_loop()
{
return COLVARS_NOT_IMPLEMENTED;
}
/// Distribute calculation of biases across threads
virtual int smp_biases_loop()
{
return COLVARS_NOT_IMPLEMENTED;
}
/// Distribute calculation of biases across threads 2nd through last, with all scripted biased on 1st thread
virtual int smp_biases_script_loop()
{
return COLVARS_NOT_IMPLEMENTED;
}
/// Index of this thread
virtual int smp_thread_id()
{
return COLVARS_NOT_IMPLEMENTED;
}
/// Number of threads sharing this address space
virtual int smp_num_threads()
{
return COLVARS_NOT_IMPLEMENTED;
}
/// Lock the proxy's shared data for access by a thread, if threads are implemented; if not implemented, does nothing
virtual int smp_lock()
{
return COLVARS_OK;
}
/// Attempt to lock the proxy's shared data
virtual int smp_trylock()
{
return COLVARS_OK;
}
/// Release the lock
virtual int smp_unlock()
{
return COLVARS_OK;
}
// **************** MULTIPLE REPLICAS COMMUNICATION ****************
// Replica exchange commands:
/// \brief Indicate if multi-replica support is available and active
virtual bool replica_enabled() { return false; }
/// \brief Index of this replica
virtual int replica_index() { return 0; }
/// \brief Total number of replica
virtual int replica_num() { return 1; }
/// \brief Synchronize replica
virtual void replica_comm_barrier() {}
/// \brief Receive data from other replica
virtual int replica_comm_recv(char* msg_data, int buf_len, int src_rep) {
return COLVARS_NOT_IMPLEMENTED;
}
/// \brief Send data to other replica
virtual int replica_comm_send(char* msg_data, int msg_len, int dest_rep) {
return COLVARS_NOT_IMPLEMENTED;
}
// **************** SCRIPTING INTERFACE ****************
/// Pointer to the scripting interface object
/// (does not need to be allocated in a new interface)
colvarscript *script;
/// is a user force script defined?
bool force_script_defined;
/// Do we have a scripting interface?
bool have_scripts;
/// Run a user-defined colvar forces script
virtual int run_force_callback() { return COLVARS_NOT_IMPLEMENTED; }
virtual int run_colvar_callback(std::string const &name,
std::vector<const colvarvalue *> const &cvcs,
colvarvalue &value)
{ return COLVARS_NOT_IMPLEMENTED; }
virtual int run_colvar_gradient_callback(std::string const &name,
std::vector<const colvarvalue *> const &cvcs,
std::vector<cvm::matrix2d<cvm::real> > &gradient)
{ return COLVARS_NOT_IMPLEMENTED; }
// **************** INPUT/OUTPUT ****************
/// Print a message to the main log
virtual void log(std::string const &message) = 0;
/// Print a message to the main log and let the rest of the program handle the error
virtual void error(std::string const &message) = 0;
/// Print a message to the main log and exit with error code
virtual void fatal_error(std::string const &message) = 0;
/// Print a message to the main log and exit normally
virtual void exit(std::string const &message)
{
cvm::error("Error: exiting without error is not implemented, returning error code.\n",
COLVARS_NOT_IMPLEMENTED);
}
// TODO the following definitions may be moved to a .cpp file
/// \brief Returns a reference to the given output channel;
/// if this is not open already, then open it
virtual std::ostream * output_stream(std::string const &output_name)
{
std::list<std::ostream *>::iterator osi = output_files.begin();
std::list<std::string>::iterator osni = output_stream_names.begin();
for ( ; osi != output_files.end(); osi++, osni++) {
if (*osni == output_name) {
return *osi;
}
}
output_stream_names.push_back(output_name);
std::ofstream * os = new std::ofstream(output_name.c_str());
if (!os->is_open()) {
cvm::error("Error: cannot write to file \""+output_name+"\".\n",
FILE_ERROR);
}
output_files.push_back(os);
return os;
}
/// \brief Closes the given output channel
virtual int close_output_stream(std::string const &output_name)
{
std::list<std::ostream *>::iterator osi = output_files.begin();
std::list<std::string>::iterator osni = output_stream_names.begin();
for ( ; osi != output_files.end(); osi++, osni++) {
if (*osni == output_name) {
((std::ofstream *) (*osi))->close();
output_files.erase(osi);
output_stream_names.erase(osni);
return COLVARS_OK;
}
}
cvm::error("Error: trying to close an output file or stream that wasn't open.\n",
BUG_ERROR);
return COLVARS_ERROR;
}
/// \brief Rename the given file, before overwriting it
virtual int backup_file(char const *filename)
{
return COLVARS_NOT_IMPLEMENTED;
}
// **************** ACCESS SYSTEM DATA ****************
/// Pass restraint energy value for current timestep to MD engine
virtual void add_energy(cvm::real energy) = 0;
/// Tell the proxy whether total forces are needed (may not always be available)
virtual void request_total_force(bool yesno)
{
if (yesno == true)
cvm::error("Error: total forces are currently not implemented.\n",
COLVARS_NOT_IMPLEMENTED);
}
/// Are total forces being used?
virtual bool total_forces_enabled() const
{
return false;
}
/// \brief Get the PBC-aware distance vector between two positions
virtual cvm::rvector position_distance(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2) = 0;
@ -346,107 +73,72 @@ public:
/// \brief Get the PBC-aware square distance between two positions;
/// may need to be reimplemented independently from position_distance() for optimization purposes
virtual cvm::real position_dist2(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2)
{
return (position_distance(pos1, pos2)).norm2();
}
cvm::atom_pos const &pos2);
/// \brief Get the closest periodic image to a reference position
/// \param pos The position to look for the closest periodic image
/// \param ref_pos The reference position
virtual void select_closest_image(cvm::atom_pos &pos,
cvm::atom_pos const &ref_pos)
{
pos = position_distance(ref_pos, pos) + ref_pos;
}
/// Tell the proxy whether total forces are needed (may not always be available)
virtual void request_total_force(bool yesno);
/// \brief Perform select_closest_image() on a set of atomic positions
///
/// After that, distance vectors can then be calculated directly,
/// without using position_distance()
void select_closest_images(std::vector<cvm::atom_pos> &pos,
cvm::atom_pos const &ref_pos)
{
for (std::vector<cvm::atom_pos>::iterator pi = pos.begin();
pi != pos.end(); ++pi) {
select_closest_image(*pi, ref_pos);
}
}
/// Are total forces being used?
virtual bool total_forces_enabled() const;
};
// **************** ACCESS ATOMIC DATA ****************
protected:
/// \brief Array of 0-based integers used to uniquely associate atoms
/// within the host program
std::vector<int> atoms_ids;
/// \brief Keep track of how many times each atom is used by a separate colvar object
std::vector<size_t> atoms_ncopies;
/// \brief Masses of the atoms (allow redefinition during a run, as done e.g. in LAMMPS)
std::vector<cvm::real> atoms_masses;
/// \brief Charges of the atoms (allow redefinition during a run, as done e.g. in LAMMPS)
std::vector<cvm::real> atoms_charges;
/// \brief Current three-dimensional positions of the atoms
std::vector<cvm::rvector> atoms_positions;
/// \brief Most recent total forces on each atom
std::vector<cvm::rvector> atoms_total_forces;
/// \brief Forces applied from colvars, to be communicated to the MD integrator
std::vector<cvm::rvector> atoms_new_colvar_forces;
/// Used by all init_atom() functions: create a slot for an atom not requested yet
inline int add_atom_slot(int atom_id)
{
atoms_ids.push_back(atom_id);
atoms_ncopies.push_back(1);
atoms_masses.push_back(1.0);
atoms_charges.push_back(0.0);
atoms_positions.push_back(cvm::rvector(0.0, 0.0, 0.0));
atoms_total_forces.push_back(cvm::rvector(0.0, 0.0, 0.0));
atoms_new_colvar_forces.push_back(cvm::rvector(0.0, 0.0, 0.0));
return (atoms_ids.size() - 1);
}
/// \brief Container of atomic data for processing by Colvars
class colvarproxy_atoms {
public:
/// Prepare this atom for collective variables calculation, selecting it by numeric index (1-based)
/// Constructor
colvarproxy_atoms();
/// Destructor
virtual ~colvarproxy_atoms();
/// Prepare this atom for collective variables calculation, selecting it by
/// numeric index (1-based)
virtual int init_atom(int atom_number) = 0;
/// Check that this atom number is valid, but do not initialize the corresponding atom yet
/// Check that this atom number is valid, but do not initialize the
/// corresponding atom yet
virtual int check_atom_id(int atom_number) = 0;
/// Select this atom for collective variables calculation, using name and residue number.
/// Not all programs support this: leave this function as is in those cases.
/// Select this atom for collective variables calculation, using name and
/// residue number. Not all programs support this: leave this function as
/// is in those cases.
virtual int init_atom(cvm::residue_id const &residue,
std::string const &atom_name,
std::string const &segment_id)
{
cvm::error("Error: initializing an atom by name and residue number is currently not supported.\n",
COLVARS_NOT_IMPLEMENTED);
return COLVARS_NOT_IMPLEMENTED;
}
std::string const &segment_id);
/// Check that this atom is valid, but do not initialize it yet
virtual int check_atom_id(cvm::residue_id const &residue,
std::string const &atom_name,
std::string const &segment_id)
{
cvm::error("Error: initializing an atom by name and residue number is currently not supported.\n",
COLVARS_NOT_IMPLEMENTED);
return COLVARS_NOT_IMPLEMENTED;
}
std::string const &segment_id);
/// \brief Used by the atom class destructor: rather than deleting the array slot
/// (costly) set the corresponding atoms_ncopies to zero
virtual void clear_atom(int index)
{
if (((size_t) index) >= atoms_ids.size()) {
cvm::error("Error: trying to disable an atom that was not previously requested.\n",
INPUT_ERROR);
}
if (atoms_ncopies[index] > 0) {
atoms_ncopies[index] -= 1;
}
}
virtual void clear_atom(int index);
/// \brief Read atom identifiers from a file \param filename name of
/// the file (usually a PDB) \param atoms array to which atoms read
/// from "filename" will be appended \param pdb_field (optiona) if
/// "filename" is a PDB file, use this field to determine which are
/// the atoms to be set
virtual int load_atoms(char const *filename,
cvm::atom_group &atoms,
std::string const &pdb_field,
double const pdb_field_value = 0.0);
/// \brief Load the coordinates for a group of atoms from a file
/// (usually a PDB); if "pos" is already allocated, the number of its
/// elements must match the number of atoms in "filename"
virtual int load_coords(char const *filename,
std::vector<cvm::atom_pos> &pos,
const std::vector<int> &indices,
std::string const &pdb_field,
double const pdb_field_value = 0.0);
/// Clear atomic data
int reset();
/// Get the numeric ID of the given atom (for the program)
inline int get_atom_id(int index) const
@ -485,120 +177,95 @@ public:
}
/// Read the current velocity of the given atom
virtual cvm::rvector get_atom_velocity(int index)
inline cvm::rvector get_atom_velocity(int index)
{
cvm::error("Error: reading the current velocity of an atom is not yet implemented.\n",
cvm::error("Error: reading the current velocity of an atom "
"is not yet implemented.\n",
COLVARS_NOT_IMPLEMENTED);
return cvm::rvector(0.0);
}
// useful functions for data management outside this class
inline std::vector<int> *modify_atom_ids() { return &atoms_ids; }
inline std::vector<cvm::real> *modify_atom_masses() { return &atoms_masses; }
inline std::vector<cvm::real> *modify_atom_charges() { return &atoms_charges; }
inline std::vector<cvm::rvector> *modify_atom_positions() { return &atoms_positions; }
inline std::vector<cvm::rvector> *modify_atom_total_forces() { return &atoms_total_forces; }
inline std::vector<cvm::rvector> *modify_atom_new_colvar_forces() { return &atoms_new_colvar_forces; }
/// \brief Read atom identifiers from a file \param filename name of
/// the file (usually a PDB) \param atoms array to which atoms read
/// from "filename" will be appended \param pdb_field (optiona) if
/// "filename" is a PDB file, use this field to determine which are
/// the atoms to be set
virtual int load_atoms(char const *filename,
cvm::atom_group &atoms,
std::string const &pdb_field,
double const pdb_field_value = 0.0)
inline std::vector<int> *modify_atom_ids()
{
cvm::error("Error: loading atom identifiers from a file is currently not implemented.\n",
COLVARS_NOT_IMPLEMENTED);
return COLVARS_NOT_IMPLEMENTED;
return &atoms_ids;
}
/// \brief Load the coordinates for a group of atoms from a file
/// (usually a PDB); if "pos" is already allocated, the number of its
/// elements must match the number of atoms in "filename"
virtual int load_coords(char const *filename,
std::vector<cvm::atom_pos> &pos,
const std::vector<int> &indices,
std::string const &pdb_field,
double const pdb_field_value = 0.0)
inline std::vector<cvm::real> *modify_atom_masses()
{
cvm::error("Error: loading atomic coordinates from a file is currently not implemented.\n");
return COLVARS_NOT_IMPLEMENTED;
return &atoms_masses;
}
// **************** ACCESS GROUP DATA ****************
inline std::vector<cvm::real> *modify_atom_charges()
{
return &atoms_charges;
}
inline std::vector<cvm::rvector> *modify_atom_positions()
{
return &atoms_positions;
}
inline std::vector<cvm::rvector> *modify_atom_total_forces()
{
return &atoms_total_forces;
}
inline std::vector<cvm::rvector> *modify_atom_new_colvar_forces()
{
return &atoms_new_colvar_forces;
}
protected:
/// \brief Array of 0-based integers used to uniquely associate atom groups
/// \brief Array of 0-based integers used to uniquely associate atoms
/// within the host program
std::vector<int> atom_groups_ids;
/// \brief Keep track of how many times each group is used by a separate cvc
std::vector<size_t> atom_groups_ncopies;
/// \brief Total masses of the atom groups
std::vector<cvm::real> atom_groups_masses;
/// \brief Total charges of the atom groups (allow redefinition during a run, as done e.g. in LAMMPS)
std::vector<cvm::real> atom_groups_charges;
/// \brief Current centers of mass of the atom groups
std::vector<cvm::rvector> atom_groups_coms;
/// \brief Most recently updated total forces on the com of each group
std::vector<cvm::rvector> atom_groups_total_forces;
std::vector<int> atoms_ids;
/// \brief Keep track of how many times each atom is used by a separate colvar object
std::vector<size_t> atoms_ncopies;
/// \brief Masses of the atoms (allow redefinition during a run, as done e.g. in LAMMPS)
std::vector<cvm::real> atoms_masses;
/// \brief Charges of the atoms (allow redefinition during a run, as done e.g. in LAMMPS)
std::vector<cvm::real> atoms_charges;
/// \brief Current three-dimensional positions of the atoms
std::vector<cvm::rvector> atoms_positions;
/// \brief Most recent total forces on each atom
std::vector<cvm::rvector> atoms_total_forces;
/// \brief Forces applied from colvars, to be communicated to the MD integrator
std::vector<cvm::rvector> atom_groups_new_colvar_forces;
std::vector<cvm::rvector> atoms_new_colvar_forces;
/// TODO Add here containers of handles to cvc objects that are computed in parallel
/// Used by all init_atom() functions: create a slot for an atom not
/// requested yet; returns the index in the arrays
int add_atom_slot(int atom_id);
};
/// \brief Container of atom group data (allow collection of aggregated atomic
/// data)
class colvarproxy_atom_groups {
public:
/// \brief Whether this proxy implementation has capability for scalable groups
virtual int scalable_group_coms()
{
return COLVARS_NOT_IMPLEMENTED;
}
/// Contructor
colvarproxy_atom_groups();
/// Used by all init_atom_group() functions: create a slot for an atom group not requested yet
// TODO Add a handle to cvc objects
inline int add_atom_group_slot(int atom_group_id)
{
atom_groups_ids.push_back(atom_group_id);
atom_groups_ncopies.push_back(1);
atom_groups_masses.push_back(1.0);
atom_groups_charges.push_back(0.0);
atom_groups_coms.push_back(cvm::rvector(0.0, 0.0, 0.0));
atom_groups_total_forces.push_back(cvm::rvector(0.0, 0.0, 0.0));
atom_groups_new_colvar_forces.push_back(cvm::rvector(0.0, 0.0, 0.0));
return (atom_groups_ids.size() - 1);
}
/// Destructor
virtual ~colvarproxy_atom_groups();
/// Clear atom group data
int reset();
/// \brief Whether this proxy implementation has capability for scalable groups
virtual int scalable_group_coms();
/// Prepare this group for collective variables calculation, selecting atoms by internal ids (0-based)
virtual int init_atom_group(std::vector<int> const &atoms_ids) // TODO Add a handle to cvc objects
{
cvm::error("Error: initializing a group outside of the colvars module is currently not supported.\n",
COLVARS_NOT_IMPLEMENTED);
return COLVARS_NOT_IMPLEMENTED;
}
virtual int init_atom_group(std::vector<int> const &atoms_ids);
/// \brief Used by the atom_group class destructor
virtual void clear_atom_group(int index)
{
if (cvm::debug()) {
log("Trying to remove/disable atom group number "+cvm::to_str(index)+"\n");
}
if (((size_t) index) >= atom_groups_ids.size()) {
cvm::error("Error: trying to disable an atom group that was not previously requested.\n",
INPUT_ERROR);
}
if (atom_groups_ncopies[index] > 0) {
atom_groups_ncopies[index] -= 1;
}
}
virtual void clear_atom_group(int index);
/// Get the numeric ID of the given atom group (for the MD program)
inline cvm::real get_atom_group_id(int index) const
inline int get_atom_group_id(int index) const
{
return atom_groups_ids[index];
}
@ -634,13 +301,288 @@ public:
}
/// Read the current velocity of the given atom group
virtual cvm::rvector get_atom_group_velocity(int index)
inline cvm::rvector get_atom_group_velocity(int index)
{
cvm::error("Error: reading the current velocity of an atom group is not yet implemented.\n",
COLVARS_NOT_IMPLEMENTED);
return cvm::rvector(0.0);
}
protected:
/// \brief Array of 0-based integers used to uniquely associate atom groups
/// within the host program
std::vector<int> atom_groups_ids;
/// \brief Keep track of how many times each group is used by a separate cvc
std::vector<size_t> atom_groups_ncopies;
/// \brief Total masses of the atom groups
std::vector<cvm::real> atom_groups_masses;
/// \brief Total charges of the atom groups (allow redefinition during a run, as done e.g. in LAMMPS)
std::vector<cvm::real> atom_groups_charges;
/// \brief Current centers of mass of the atom groups
std::vector<cvm::rvector> atom_groups_coms;
/// \brief Most recently updated total forces on the com of each group
std::vector<cvm::rvector> atom_groups_total_forces;
/// \brief Forces applied from colvars, to be communicated to the MD integrator
std::vector<cvm::rvector> atom_groups_new_colvar_forces;
/// Used by all init_atom_group() functions: create a slot for an atom group not requested yet
int add_atom_group_slot(int atom_group_id);
};
/// \brief Methods for SMP parallelization
class colvarproxy_smp {
public:
/// Constructor
colvarproxy_smp();
/// Destructor
virtual ~colvarproxy_smp();
/// Whether threaded parallelization should be used (TODO: make this a
/// cvm::deps feature)
bool b_smp_active;
/// Whether threaded parallelization is available (TODO: make this a cvm::deps feature)
virtual int smp_enabled();
/// Distribute calculation of colvars (and their components) across threads
virtual int smp_colvars_loop();
/// Distribute calculation of biases across threads
virtual int smp_biases_loop();
/// Distribute calculation of biases across threads 2nd through last, with all scripted biased on 1st thread
virtual int smp_biases_script_loop();
/// Index of this thread
virtual int smp_thread_id();
/// Number of threads sharing this address space
virtual int smp_num_threads();
/// Lock the proxy's shared data for access by a thread, if threads are implemented; if not implemented, does nothing
virtual int smp_lock();
/// Attempt to lock the proxy's shared data
virtual int smp_trylock();
/// Release the lock
virtual int smp_unlock();
};
/// \brief Methods for multiple-replica communication
class colvarproxy_replicas {
public:
/// Constructor
colvarproxy_replicas();
/// Destructor
virtual ~colvarproxy_replicas();
/// \brief Indicate if multi-replica support is available and active
virtual bool replica_enabled();
/// \brief Index of this replica
virtual int replica_index();
/// \brief Total number of replica
virtual int replica_num();
/// \brief Synchronize replica
virtual void replica_comm_barrier();
/// \brief Receive data from other replica
virtual int replica_comm_recv(char* msg_data, int buf_len, int src_rep);
/// \brief Send data to other replica
virtual int replica_comm_send(char* msg_data, int msg_len, int dest_rep);
};
/// Method for scripting language interface (Tcl or Python)
class colvarproxy_script {
public:
/// Constructor
colvarproxy_script();
/// Destructor
virtual ~colvarproxy_script();
/// Convert a script object (Tcl or Python call argument) to a C string
virtual char *script_obj_to_str(unsigned char *obj);
/// Pointer to the scripting interface object
/// (does not need to be allocated in a new interface)
colvarscript *script;
/// is a user force script defined?
bool force_script_defined;
/// Do we have a scripting interface?
bool have_scripts;
/// Run a user-defined colvar forces script
virtual int run_force_callback();
virtual int run_colvar_callback(
std::string const &name,
std::vector<const colvarvalue *> const &cvcs,
colvarvalue &value);
virtual int run_colvar_gradient_callback(
std::string const &name,
std::vector<const colvarvalue *> const &cvcs,
std::vector<cvm::matrix2d<cvm::real> > &gradient);
};
/// Methods for data input/output
class colvarproxy_io {
public:
/// Constructor
colvarproxy_io();
/// Destructor
virtual ~colvarproxy_io();
/// \brief Save the current frame number in the argument given
// Returns error code
virtual int get_frame(long int &);
/// \brief Set the current frame number (as well as colvarmodule::it)
// Returns error code
virtual int set_frame(long int);
/// \brief Returns a reference to the given output channel;
/// if this is not open already, then open it
virtual std::ostream *output_stream(std::string const &output_name,
std::ios_base::openmode mode =
std::ios_base::out);
/// \brief Flushes the given output channel
virtual int flush_output_stream(std::ostream *os);
/// \brief Closes the given output channel
virtual int close_output_stream(std::string const &output_name);
/// \brief Rename the given file, before overwriting it
virtual int backup_file(char const *filename);
/// \brief Rename the given file, before overwriting it
inline int backup_file(std::string const &filename)
{
return backup_file(filename.c_str());
}
/// \brief Prefix of the input state file
inline std::string & input_prefix()
{
return input_prefix_str;
}
/// \brief Prefix to be used for output restart files
inline std::string & restart_output_prefix()
{
return restart_output_prefix_str;
}
/// \brief Prefix to be used for output files (final system
/// configuration)
inline std::string & output_prefix()
{
return output_prefix_str;
}
protected:
/// \brief Prefix to be used for input files (restarts, not
/// configuration)
std::string input_prefix_str, output_prefix_str, restart_output_prefix_str;
/// \brief Currently opened output files: by default, these are ofstream objects.
/// Allows redefinition to implement different output mechanisms
std::list<std::ostream *> output_files;
/// \brief Identifiers for output_stream objects: by default, these are the names of the files
std::list<std::string> output_stream_names;
};
/// \brief Interface between the collective variables module and
/// the simulation or analysis program (NAMD, VMD, LAMMPS...).
/// This is the base class: each interfaced program is supported by a derived class.
/// Only pure virtual functions ("= 0") must be reimplemented to ensure baseline functionality.
class colvarproxy
: public colvarproxy_system,
public colvarproxy_atoms,
public colvarproxy_atom_groups,
public colvarproxy_smp,
public colvarproxy_replicas,
public colvarproxy_script,
public colvarproxy_io
{
public:
/// Pointer to the main object
colvarmodule *colvars;
/// Constructor
colvarproxy();
/// Destructor
virtual ~colvarproxy();
/// \brief Reset proxy state, e.g. requested atoms
virtual int reset();
/// (Re)initialize required member data after construction
virtual int setup();
/// \brief Update data required by the colvars module (e.g. cache atom positions)
///
/// TODO Break up colvarproxy_namd and colvarproxy_lammps function into these
virtual int update_input();
/// \brief Update data based from the results of a module update (e.g. send forces)
virtual int update_output();
/// Print a message to the main log
virtual void log(std::string const &message) = 0;
/// Print a message to the main log and let the rest of the program handle the error
virtual void error(std::string const &message) = 0;
/// Print a message to the main log and exit with error code
virtual void fatal_error(std::string const &message) = 0;
/// \brief Restarts will be written each time this number of steps has passed
virtual size_t restart_frequency();
/// Whether a simulation is running (warn against irrecovarable errors)
inline bool simulation_running() const
{
return b_simulation_running;
}
protected:
/// Whether a simulation is running (warn against irrecovarable errors)
bool b_simulation_running;
};

8
lib/colvars/colvars_version.h Normal file
View File

@ -0,0 +1,8 @@
#define COLVARS_VERSION "2017-07-15"
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/colvars/colvars
// Please update all Colvars source files before making any changes.
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.

150
lib/colvars/colvarscript.cpp
View File

@ -12,6 +12,7 @@
#include <string.h>
#include "colvarscript.h"
#include "colvarproxy.h"
#include "colvardeps.h"
@ -27,7 +28,7 @@ extern "C" {
// Generic hooks; NAMD and VMD have Tcl-specific versions in the respective proxies
int run_colvarscript_command(int argc, const char **argv)
int run_colvarscript_command(int objc, unsigned char *const objv[])
{
colvarproxy *cvp = cvm::proxy;
if (!cvp) {
@ -37,7 +38,7 @@ extern "C" {
cvm::error("Called run_colvarscript_command without a script object initialized.\n");
return -1;
}
return cvp->script->run(argc, argv);
return cvp->script->run(objc, objv);
}
const char * get_colvarscript_result()
@ -53,30 +54,52 @@ extern "C" {
/// Run method based on given arguments
int colvarscript::run(int argc, char const *argv[]) {
result = "";
int colvarscript::run(int objc, unsigned char *const objv[])
{
result.clear();
if (cvm::debug()) {
cvm::log("Called script run with " + cvm::to_str(argc) + " args");
for (int i = 0; i < argc; i++) { cvm::log(argv[i]); }
cvm::log("Called script run with " + cvm::to_str(objc) + " args:");
for (int i = 0; i < objc; i++) {
cvm::log(obj_to_str(objv[i]));
}
}
if (argc < 2) {
if (objc < 2) {
result = help_string();
return COLVARS_OK;
}
std::string cmd = argv[1];
std::string const cmd(obj_to_str(objv[1]));
int error_code = COLVARS_OK;
if (cmd == "colvar") {
return proc_colvar(argc-1, &(argv[1]));
if (objc < 3) {
result = "Missing parameters\n" + help_string();
return COLVARSCRIPT_ERROR;
}
std::string const name(obj_to_str(objv[2]));
colvar *cv = cvm::colvar_by_name(name);
if (cv == NULL) {
result = "Colvar not found: " + name;
return COLVARSCRIPT_ERROR;
}
return proc_colvar(cv, objc-1, &(objv[1]));
}
if (cmd == "bias") {
return proc_bias(argc-1, &(argv[1]));
if (objc < 3) {
result = "Missing parameters\n" + help_string();
return COLVARSCRIPT_ERROR;
}
std::string const name(obj_to_str(objv[2]));
colvarbias *b = cvm::bias_by_name(name);
if (b == NULL) {
result = "Bias not found: " + name;
return COLVARSCRIPT_ERROR;
}
return proc_bias(b, objc-1, &(objv[1]));
}
if (cmd == "version") {
@ -102,20 +125,20 @@ int colvarscript::run(int argc, char const *argv[]) {
error_code |= colvars->calc();
error_code |= proxy->update_output();
if (error_code) {
result += "Error updating the colvars module.\n";
result += "Error updating the Colvars module.\n";
}
return error_code;
}
if (cmd == "list") {
if (argc == 2) {
if (objc == 2) {
for (std::vector<colvar *>::iterator cvi = colvars->colvars.begin();
cvi != colvars->colvars.end();
++cvi) {
result += (cvi == colvars->colvars.begin() ? "" : " ") + (*cvi)->name;
}
return COLVARS_OK;
} else if (argc == 3 && !strcmp(argv[2], "biases")) {
} else if (objc == 3 && !strcmp(obj_to_str(objv[2]), "biases")) {
for (std::vector<colvarbias *>::iterator bi = colvars->biases.begin();
bi != colvars->biases.end();
++bi) {
@ -130,11 +153,11 @@ int colvarscript::run(int argc, char const *argv[]) {
/// Parse config from file
if (cmd == "configfile") {
if (argc < 3) {
if (objc < 3) {
result = "Missing arguments\n" + help_string();
return COLVARSCRIPT_ERROR;
}
if (colvars->read_config_file(argv[2]) == COLVARS_OK) {
if (colvars->read_config_file(obj_to_str(objv[2])) == COLVARS_OK) {
return COLVARS_OK;
} else {
result = "Error parsing configuration file";
@ -144,11 +167,11 @@ int colvarscript::run(int argc, char const *argv[]) {
/// Parse config from string
if (cmd == "config") {
if (argc < 3) {
if (objc < 3) {
result = "Missing arguments\n" + help_string();
return COLVARSCRIPT_ERROR;
}
std::string conf = argv[2];
std::string const conf(obj_to_str(objv[2]));
if (colvars->read_config_string(conf) == COLVARS_OK) {
return COLVARS_OK;
} else {
@ -159,11 +182,11 @@ int colvarscript::run(int argc, char const *argv[]) {
/// Load an input state file
if (cmd == "load") {
if (argc < 3) {
if (objc < 3) {
result = "Missing arguments\n" + help_string();
return COLVARSCRIPT_ERROR;
}
proxy->input_prefix() = argv[2];
proxy->input_prefix() = obj_to_str(objv[2]);
if (colvars->setup_input() == COLVARS_OK) {
return COLVARS_OK;
} else {
@ -174,11 +197,11 @@ int colvarscript::run(int argc, char const *argv[]) {
/// Save to an output state file
if (cmd == "save") {
if (argc < 3) {
if (objc < 3) {
result = "Missing arguments";
return COLVARSCRIPT_ERROR;
}
proxy->output_prefix_str = argv[2];
proxy->output_prefix() = obj_to_str(objv[2]);
int error = 0;
error |= colvars->setup_output();
error |= colvars->write_output_files();
@ -200,7 +223,7 @@ int colvarscript::run(int argc, char const *argv[]) {
}
if (cmd == "frame") {
if (argc == 2) {
if (objc == 2) {
long int f;
int error = proxy->get_frame(f);
if (error == COLVARS_OK) {
@ -210,10 +233,10 @@ int colvarscript::run(int argc, char const *argv[]) {
result = "Frame number is not available";
return COLVARSCRIPT_ERROR;
}
} else if (argc == 3) {
} else if (objc == 3) {
// Failure of this function does not trigger an error, but
// returns nonzero, to let scripts detect available frames
int error = proxy->set_frame(strtol(argv[2], NULL, 10));
int error = proxy->set_frame(strtol(obj_to_str(objv[2]), NULL, 10));
result = cvm::to_str(error == COLVARS_OK ? 0 : -1);
return COLVARS_OK;
} else {
@ -223,8 +246,8 @@ int colvarscript::run(int argc, char const *argv[]) {
}
if (cmd == "addenergy") {
if (argc == 3) {
colvars->total_bias_energy += strtod(argv[2], NULL);
if (objc == 3) {
colvars->total_bias_energy += strtod(obj_to_str(objv[2]), NULL);
return COLVARS_OK;
} else {
result = "Wrong arguments to command \"addenergy\"\n" + help_string();
@ -237,19 +260,9 @@ int colvarscript::run(int argc, char const *argv[]) {
}
int colvarscript::proc_colvar(int argc, char const *argv[]) {
if (argc < 3) {
result = "Missing parameters\n" + help_string();
return COLVARSCRIPT_ERROR;
}
int colvarscript::proc_colvar(colvar *cv, int objc, unsigned char *const objv[]) {
std::string name = argv[1];
colvar *cv = cvm::colvar_by_name(name);
if (cv == NULL) {
result = "Colvar not found: " + name;
return COLVARSCRIPT_ERROR;
}
std::string subcmd = argv[2];
std::string const subcmd(obj_to_str(objv[2]));
if (subcmd == "value") {
result = (cv->value()).to_simple_string();
@ -278,11 +291,11 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
for (i = 0; i < cv->biases.size(); i++) {
delete cv->biases[i];
}
cv->biases.resize(0);
cv->biases.clear();
// colvar destructor is tasked with the cleanup
delete cv;
// TODO this could be done by the destructors
colvars->write_traj_label(colvars->cv_traj_os);
colvars->write_traj_label(*(colvars->cv_traj_os));
return COLVARS_OK;
}
@ -308,11 +321,11 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
}
if (subcmd == "addforce") {
if (argc < 4) {
if (objc < 4) {
result = "addforce: missing parameter: force value\n" + help_string();
return COLVARSCRIPT_ERROR;
}
std::string f_str = argv[3];
std::string const f_str(obj_to_str(objv[3]));
std::istringstream is(f_str);
is.width(cvm::cv_width);
is.precision(cvm::cv_prec);
@ -328,11 +341,11 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
}
if (subcmd == "cvcflags") {
if (argc < 4) {
if (objc < 4) {
result = "cvcflags: missing parameter: vector of flags";
return COLVARSCRIPT_ERROR;
}
std::string flags_str = argv[3];
std::string const flags_str(obj_to_str(objv[3]));
std::istringstream is(flags_str);
std::vector<bool> flags;
@ -351,7 +364,7 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
}
if ((subcmd == "get") || (subcmd == "set") || (subcmd == "state")) {
return proc_features(cv, argc, argv);
return proc_features(cv, objc, objv);
}
result = "Syntax error\n" + help_string();
@ -359,20 +372,10 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
}
int colvarscript::proc_bias(int argc, char const *argv[]) {
if (argc < 3) {
result = "Missing parameters\n" + help_string();
return COLVARSCRIPT_ERROR;
}
int colvarscript::proc_bias(colvarbias *b, int objc, unsigned char *const objv[]) {
std::string name = argv[1];
colvarbias *b = cvm::bias_by_name(name);
if (b == NULL) {
result = "Bias not found: " + name;
return COLVARSCRIPT_ERROR;
}
std::string subcmd = argv[2];
std::string const key(obj_to_str(objv[0]));
std::string const subcmd(obj_to_str(objv[2]));
if (subcmd == "energy") {
result = cvm::to_str(b->get_energy());
@ -422,16 +425,16 @@ int colvarscript::proc_bias(int argc, char const *argv[]) {
// the bias destructor takes care of the cleanup at cvm level
delete b;
// TODO this could be done by the destructors
colvars->write_traj_label(colvars->cv_traj_os);
colvars->write_traj_label(*(colvars->cv_traj_os));
return COLVARS_OK;
}
if ((subcmd == "get") || (subcmd == "set") || (subcmd == "state")) {
return proc_features(b, argc, argv);
return proc_features(b, objc, objv);
}
if (argc >= 4) {
std::string param = argv[3];
if (objc >= 4) {
std::string const param(obj_to_str(objv[3]));
if (subcmd == "count") {
int index;
if (!(std::istringstream(param) >> index)) {
@ -452,11 +455,11 @@ int colvarscript::proc_bias(int argc, char const *argv[]) {
int colvarscript::proc_features(colvardeps *obj,
int argc, char const *argv[]) {
int objc, unsigned char *const objv[]) {
// size was already checked before calling
std::string subcmd = argv[2];
std::string const subcmd(obj_to_str(objv[2]));
if (argc == 3) {
if (objc == 3) {
if (subcmd == "state") {
// TODO make this returned as result?
obj->print_state();
@ -470,7 +473,7 @@ int colvarscript::proc_features(colvardeps *obj,
if ((subcmd == "get") || (subcmd == "set")) {
std::vector<colvardeps::feature *> &features = obj->features();
std::string const req_feature(argv[3]);
std::string const req_feature(obj_to_str(objv[3]));
colvardeps::feature *f = NULL;
int fid = 0;
for (fid = 0; fid < int(features.size()); fid++) {
@ -499,9 +502,9 @@ int colvarscript::proc_features(colvardeps *obj,
}
if (subcmd == "set") {
if (argc == 5) {
if (objc == 5) {
std::string const yesno =
colvarparse::to_lower_cppstr(std::string(argv[4]));
colvarparse::to_lower_cppstr(std::string(obj_to_str(objv[4])));
if ((yesno == std::string("yes")) ||
(yesno == std::string("on")) ||
(yesno == std::string("1"))) {
@ -510,10 +513,7 @@ int colvarscript::proc_features(colvardeps *obj,
} else if ((yesno == std::string("no")) ||
(yesno == std::string("off")) ||
(yesno == std::string("0"))) {
// TODO disable() function does not exist yet,
// dependencies will not be resolved
// obj->disable(fid);
obj->set_enabled(fid, false);
obj->disable(fid);
return COLVARS_OK;
}
}
@ -533,11 +533,11 @@ std::string colvarscript::help_string()
std::string buf;
buf = "Usage: cv <subcommand> [args...]\n\
\n\
Managing the colvars module:\n\
Managing the Colvars module:\n\
configfile <file name> -- read configuration from a file\n\
config <string> -- read configuration from the given string\n\
reset -- delete all internal configuration\n\
delete -- delete this colvars module instance\n\
delete -- delete this Colvars module instance\n\
version -- return version of colvars code\n\
\n\
Input and output:\n\

20
lib/colvars/colvarscript.h
View File

@ -41,22 +41,30 @@ public:
/// If an error is returned by one of the methods, it should set this to the error message
std::string result;
/// Run script command with given positional arguments
int run(int argc, char const *argv[]);
/// Run script command with given positional arguments (objects)
int run(int objc, unsigned char *const objv[]);
private:
/// Run subcommands on colvar
int proc_colvar(int argc, char const *argv[]);
int proc_colvar(colvar *cv, int argc, unsigned char *const argv[]);
/// Run subcommands on bias
int proc_bias(int argc, char const *argv[]);
int proc_bias(colvarbias *b, int argc, unsigned char *const argv[]);
/// Run subcommands on base colvardeps object (colvar, bias, ...)
int proc_features(colvardeps *obj,
int argc, char const *argv[]);
int argc, unsigned char *const argv[]);
/// Builds and return a short help
/// Build and return a short help
std::string help_string(void);
public:
inline char const *obj_to_str(unsigned char *const obj)
{
return cvm::proxy->script_obj_to_str(obj);
}
};

5
lib/colvars/colvartypes.h
View File

@ -91,6 +91,11 @@ public:
data.resize(n);
}
inline void clear()
{
data.clear();
}
inline T & operator [] (size_t const i) {
return data[i];
}

275
lib/colvars/colvarvalue.cpp
View File

@ -16,6 +16,274 @@
std::string const colvarvalue::type_desc(Type t)
{
switch (t) {
case colvarvalue::type_scalar:
return "scalar number"; break;
case colvarvalue::type_3vector:
return "3-dimensional vector"; break;
case colvarvalue::type_unit3vector:
return "3-dimensional unit vector"; break;
case colvarvalue::type_unit3vectorderiv:
return "derivative of a 3-dimensional unit vector"; break;
case colvarvalue::type_quaternion:
return "4-dimensional unit quaternion"; break;
case colvarvalue::type_quaternionderiv:
return "4-dimensional tangent vector"; break;
case colvarvalue::type_vector:
return "n-dimensional vector"; break;
case colvarvalue::type_notset:
// fallthrough
default:
return "not set"; break;
}
}
std::string const colvarvalue::type_keyword(Type t)
{
switch (t) {
case colvarvalue::type_notset:
default:
return "not_set"; break;
case colvarvalue::type_scalar:
return "scalar"; break;
case colvarvalue::type_3vector:
return "vector3"; break;
case colvarvalue::type_unit3vector:
return "unit_vector3"; break;
case colvarvalue::type_unit3vectorderiv:
return ""; break;
case colvarvalue::type_quaternion:
return "unit_quaternion"; break;
case colvarvalue::type_quaternionderiv:
return ""; break;
case colvarvalue::type_vector:
return "vector"; break;
}
}
size_t colvarvalue::num_df(Type t)
{
switch (t) {
case colvarvalue::type_notset:
default:
return 0; break;
case colvarvalue::type_scalar:
return 1; break;
case colvarvalue::type_3vector:
return 3; break;
case colvarvalue::type_unit3vector:
case colvarvalue::type_unit3vectorderiv:
return 2; break;
case colvarvalue::type_quaternion:
case colvarvalue::type_quaternionderiv:
return 3; break;
case colvarvalue::type_vector:
// the size of a vector is unknown without its object
return 0; break;
}
}
size_t colvarvalue::num_dimensions(Type t)
{
switch (t) {
case colvarvalue::type_notset:
default:
return 0; break;
case colvarvalue::type_scalar:
return 1; break;
case colvarvalue::type_3vector:
case colvarvalue::type_unit3vector:
case colvarvalue::type_unit3vectorderiv:
return 3; break;
case colvarvalue::type_quaternion:
case colvarvalue::type_quaternionderiv:
return 4; break;
case colvarvalue::type_vector:
// the size of a vector is unknown without its object
return 0; break;
}
}
void colvarvalue::reset()
{
switch (value_type) {
case colvarvalue::type_scalar:
real_value = 0.0;
break;
case colvarvalue::type_3vector:
case colvarvalue::type_unit3vector:
case colvarvalue::type_unit3vectorderiv:
rvector_value.reset();
break;
case colvarvalue::type_quaternion:
case colvarvalue::type_quaternionderiv:
quaternion_value.reset();
break;
case colvarvalue::type_vector:
vector1d_value.reset();
break;
case colvarvalue::type_notset:
default:
break;
}
}
void colvarvalue::apply_constraints()
{
switch (value_type) {
case colvarvalue::type_scalar:
case colvarvalue::type_3vector:
case colvarvalue::type_unit3vectorderiv:
case colvarvalue::type_quaternionderiv:
break;
case colvarvalue::type_unit3vector:
rvector_value /= std::sqrt(rvector_value.norm2());
break;
case colvarvalue::type_quaternion:
quaternion_value /= std::sqrt(quaternion_value.norm2());
break;
case colvarvalue::type_vector:
if (elem_types.size() > 0) {
// if we have information about non-scalar types, use it
size_t i;
for (i = 0; i < elem_types.size(); i++) {
if (elem_sizes[i] == 1) continue; // TODO this can be optimized further
colvarvalue cvtmp(vector1d_value.slice(elem_indices[i],
elem_indices[i] + elem_sizes[i]), elem_types[i]);
cvtmp.apply_constraints();
set_elem(i, cvtmp);
}
}
break;
case colvarvalue::type_notset:
default:
break;
}
}
void colvarvalue::type(Type const &vti)
{
if (vti != value_type) {
// reset the value based on the previous type
reset();
if ((value_type == type_vector) && (vti != type_vector)) {
vector1d_value.clear();
}
value_type = vti;
}
}
void colvarvalue::type(colvarvalue const &x)
{
if (x.type() != value_type) {
// reset the value based on the previous type
reset();
if (value_type == type_vector) {
vector1d_value.clear();
}
value_type = x.type();
}
if (x.type() == type_vector) {
vector1d_value.resize(x.vector1d_value.size());
}
}
void colvarvalue::is_derivative()
{
switch (value_type) {
case colvarvalue::type_scalar:
case colvarvalue::type_3vector:
case colvarvalue::type_unit3vectorderiv:
case colvarvalue::type_quaternionderiv:
break;
case colvarvalue::type_unit3vector:
type(colvarvalue::type_unit3vectorderiv);
break;
case colvarvalue::type_quaternion:
type(colvarvalue::type_quaternionderiv);
break;
case colvarvalue::type_vector:
// TODO
break;
case colvarvalue::type_notset:
default:
break;
}
}
colvarvalue::colvarvalue(colvarvalue const &x)
: value_type(x.type())
{
switch (x.type()) {
case type_scalar:
real_value = x.real_value;
break;
case type_3vector:
case type_unit3vector:
case type_unit3vectorderiv:
rvector_value = x.rvector_value;
break;
case type_quaternion:
case type_quaternionderiv:
quaternion_value = x.quaternion_value;
break;
case type_vector:
vector1d_value = x.vector1d_value;
elem_types = x.elem_types;
elem_indices = x.elem_indices;
elem_sizes = x.elem_sizes;
case type_notset:
default:
break;
}
}
colvarvalue::colvarvalue(cvm::vector1d<cvm::real> const &v, Type vti)
{
if ((vti != type_vector) && (v.size() != num_dimensions(vti))) {
cvm::error("Error: trying to initialize a variable of type \""+type_desc(vti)+
"\" using a vector of size "+cvm::to_str(v.size())+
".\n");
value_type = type_notset;
} else {
value_type = vti;
switch (vti) {
case type_scalar:
real_value = v[0];
break;
case type_3vector:
case type_unit3vector:
case type_unit3vectorderiv:
rvector_value = cvm::rvector(v);
break;
case type_quaternion:
case type_quaternionderiv:
quaternion_value = cvm::quaternion(v);
break;
case type_vector:
vector1d_value = v;
break;
case type_notset:
default:
break;
}
}
}
void colvarvalue::add_elem(colvarvalue const &x)
{
if (this->value_type != type_vector) {
@ -111,6 +379,13 @@ void colvarvalue::set_random()
}
void colvarvalue::undef_op() const
{
cvm::error("Error: Undefined operation on a colvar of type \""+
type_desc(this->type())+"\".\n");
}
// binary operations between two colvarvalues
colvarvalue operator + (colvarvalue const &x1,

322
lib/colvars/colvarvalue.h
View File

@ -169,38 +169,14 @@ public:
}
/// Set the type explicitly
inline void type(Type const &vti)
{
if (vti != value_type) {
// reset the value based on the previous type
reset();
if ((value_type == type_vector) && (vti != type_vector)) {
vector1d_value.resize(0);
}
value_type = vti;
}
}
void type(Type const &vti);
/// Set the type after another \link colvarvalue \endlink
inline void type(colvarvalue const &x)
{
if (x.type() != value_type) {
// reset the value based on the previous type
reset();
if (value_type == type_vector) {
vector1d_value.resize(0);
}
value_type = x.type();
}
if (x.type() == type_vector) {
vector1d_value.resize(x.vector1d_value.size());
}
}
void type(colvarvalue const &x);
/// Make the type a derivative of the original type
/// (so that its constraints do not apply)
inline void is_derivative();
void is_derivative();
/// Square norm of this colvarvalue
cvm::real norm2() const;
@ -303,28 +279,10 @@ public:
void set_elem(int const icv, colvarvalue const &x);
/// Get a scalar number out of an element of the vector
inline cvm::real operator [] (int const i) const
{
if (vector1d_value.size() > 0) {
return vector1d_value[i];
} else {
cvm::error("Error: trying to use as a vector a variable that is not initialized as such.\n");
return 0.0;
}
}
cvm::real operator [] (int const i) const;
/// Use an element of the vector as a scalar number
inline cvm::real & operator [] (int const i)
{
if (vector1d_value.size() > 0) {
return vector1d_value[i];
} else {
cvm::error("Error: trying to use as a vector a variable that is not initialized as such.\n");
real_value = 0.0;
return real_value;
}
}
cvm::real & operator [] (int const i);
/// Ensure that the two types are the same within a binary operator
int static check_types(colvarvalue const &x1, colvarvalue const &x2);
@ -389,101 +347,6 @@ public:
};
inline std::string const colvarvalue::type_desc(Type t)
{
switch (t) {
case colvarvalue::type_scalar:
return "scalar number"; break;
case colvarvalue::type_3vector:
return "3-dimensional vector"; break;
case colvarvalue::type_unit3vector:
return "3-dimensional unit vector"; break;
case colvarvalue::type_unit3vectorderiv:
return "derivative of a 3-dimensional unit vector"; break;
case colvarvalue::type_quaternion:
return "4-dimensional unit quaternion"; break;
case colvarvalue::type_quaternionderiv:
return "4-dimensional tangent vector"; break;
case colvarvalue::type_vector:
return "n-dimensional vector"; break;
case colvarvalue::type_notset:
// fallthrough
default:
return "not set"; break;
}
}
inline std::string const colvarvalue::type_keyword(Type t)
{
switch (t) {
case colvarvalue::type_notset:
default:
return "not_set"; break;
case colvarvalue::type_scalar:
return "scalar"; break;
case colvarvalue::type_3vector:
return "vector3"; break;
case colvarvalue::type_unit3vector:
return "unit_vector3"; break;
case colvarvalue::type_unit3vectorderiv:
return ""; break;
case colvarvalue::type_quaternion:
return "unit_quaternion"; break;
case colvarvalue::type_quaternionderiv:
return ""; break;
case colvarvalue::type_vector:
return "vector"; break;
}
}
inline size_t colvarvalue::num_df(Type t)
{
switch (t) {
case colvarvalue::type_notset:
default:
return 0; break;
case colvarvalue::type_scalar:
return 1; break;
case colvarvalue::type_3vector:
return 3; break;
case colvarvalue::type_unit3vector:
case colvarvalue::type_unit3vectorderiv:
return 2; break;
case colvarvalue::type_quaternion:
case colvarvalue::type_quaternionderiv:
return 3; break;
case colvarvalue::type_vector:
// the size of a vector is unknown without its object
return 0; break;
}
}
inline size_t colvarvalue::num_dimensions(Type t)
{
switch (t) {
case colvarvalue::type_notset:
default:
return 0; break;
case colvarvalue::type_scalar:
return 1; break;
case colvarvalue::type_3vector:
case colvarvalue::type_unit3vector:
case colvarvalue::type_unit3vectorderiv:
return 3; break;
case colvarvalue::type_quaternion:
case colvarvalue::type_quaternionderiv:
return 4; break;
case colvarvalue::type_vector:
// the size of a vector is unknown without its object
return 0; break;
}
}
inline size_t colvarvalue::size() const
{
switch (value_type) {
@ -505,62 +368,48 @@ inline size_t colvarvalue::size() const
}
inline colvarvalue::colvarvalue(colvarvalue const &x)
: value_type(x.type())
inline cvm::real colvarvalue::operator [] (int const i) const
{
switch (x.type()) {
case type_scalar:
real_value = x.real_value;
break;
case type_3vector:
case type_unit3vector:
case type_unit3vectorderiv:
rvector_value = x.rvector_value;
break;
case type_quaternion:
case type_quaternionderiv:
quaternion_value = x.quaternion_value;
break;
case type_vector:
vector1d_value = x.vector1d_value;
elem_types = x.elem_types;
elem_indices = x.elem_indices;
elem_sizes = x.elem_sizes;
case type_notset:
switch (value_type) {
case colvarvalue::type_notset:
default:
break;
cvm::error("Error: trying to access a colvar value "
"that is not initialized.\n", BUG_ERROR);
return 0.0; break;
case colvarvalue::type_scalar:
return real_value; break;
case colvarvalue::type_3vector:
case colvarvalue::type_unit3vector:
case colvarvalue::type_unit3vectorderiv:
return rvector_value[i]; break;
case colvarvalue::type_quaternion:
case colvarvalue::type_quaternionderiv:
return quaternion_value[i]; break;
case colvarvalue::type_vector:
return vector1d_value[i]; break;
}
}
inline colvarvalue::colvarvalue(cvm::vector1d<cvm::real> const &v, Type vti)
inline cvm::real & colvarvalue::operator [] (int const i)
{
if ((vti != type_vector) && (v.size() != num_dimensions(vti))) {
cvm::error("Error: trying to initialize a variable of type \""+type_desc(vti)+
"\" using a vector of size "+cvm::to_str(v.size())+
".\n");
value_type = type_notset;
} else {
value_type = vti;
switch (vti) {
case type_scalar:
real_value = v[0];
break;
case type_3vector:
case type_unit3vector:
case type_unit3vectorderiv:
rvector_value = cvm::rvector(v);
break;
case type_quaternion:
case type_quaternionderiv:
quaternion_value = cvm::quaternion(v);
break;
case type_vector:
vector1d_value = v;
break;
case type_notset:
switch (value_type) {
case colvarvalue::type_notset:
default:
break;
}
cvm::error("Error: trying to access a colvar value "
"that is not initialized.\n", BUG_ERROR);
return real_value; break;
case colvarvalue::type_scalar:
return real_value; break;
case colvarvalue::type_3vector:
case colvarvalue::type_unit3vector:
case colvarvalue::type_unit3vectorderiv:
return rvector_value[i]; break;
case colvarvalue::type_quaternion:
case colvarvalue::type_quaternionderiv:
return quaternion_value[i]; break;
case colvarvalue::type_vector:
return vector1d_value[i]; break;
}
}
@ -638,13 +487,6 @@ inline int colvarvalue::check_types_assign(colvarvalue::Type const &vt1,
}
inline void colvarvalue::undef_op() const
{
cvm::error("Error: Undefined operation on a colvar of type \""+
type_desc(this->type())+"\".\n");
}
inline colvarvalue & colvarvalue::operator = (colvarvalue const &x)
{
check_types_assign(this->type(), x.type());
@ -704,6 +546,7 @@ inline void colvarvalue::operator += (colvarvalue const &x)
}
}
inline void colvarvalue::operator -= (colvarvalue const &x)
{
colvarvalue::check_types(*this, x);
@ -802,89 +645,6 @@ inline cvm::vector1d<cvm::real> const colvarvalue::as_vector() const
}
inline void colvarvalue::reset()
{
switch (value_type) {
case colvarvalue::type_scalar:
real_value = 0.0;
break;
case colvarvalue::type_3vector:
case colvarvalue::type_unit3vector:
case colvarvalue::type_unit3vectorderiv:
rvector_value.reset();
break;
case colvarvalue::type_quaternion:
case colvarvalue::type_quaternionderiv:
quaternion_value.reset();
break;
case colvarvalue::type_vector:
vector1d_value.reset();
break;
case colvarvalue::type_notset:
default:
break;
}
}
inline void colvarvalue::apply_constraints()
{
switch (value_type) {
case colvarvalue::type_scalar:
case colvarvalue::type_3vector:
case colvarvalue::type_unit3vectorderiv:
case colvarvalue::type_quaternionderiv:
break;
case colvarvalue::type_unit3vector:
rvector_value /= std::sqrt(rvector_value.norm2());
break;
case colvarvalue::type_quaternion:
quaternion_value /= std::sqrt(quaternion_value.norm2());
break;
case colvarvalue::type_vector:
if (elem_types.size() > 0) {
// if we have information about non-scalar types, use it
size_t i;
for (i = 0; i < elem_types.size(); i++) {
if (elem_sizes[i] == 1) continue; // TODO this can be optimized further
colvarvalue cvtmp(vector1d_value.slice(elem_indices[i],
elem_indices[i] + elem_sizes[i]), elem_types[i]);
cvtmp.apply_constraints();
set_elem(i, cvtmp);
}
}
break;
case colvarvalue::type_notset:
default:
break;
}
}
inline void colvarvalue::is_derivative()
{
switch (value_type) {
case colvarvalue::type_scalar:
case colvarvalue::type_3vector:
case colvarvalue::type_unit3vectorderiv:
case colvarvalue::type_quaternionderiv:
break;
case colvarvalue::type_unit3vector:
type(colvarvalue::type_unit3vectorderiv);
break;
case colvarvalue::type_quaternion:
type(colvarvalue::type_quaternionderiv);
break;
case colvarvalue::type_vector:
// TODO
break;
case colvarvalue::type_notset:
default:
break;
}
}
inline cvm::real colvarvalue::norm2() const
{
switch (value_type) {

4
src/USER-COLVARS/colvarproxy_lammps.cpp
View File

@ -308,10 +308,6 @@ void colvarproxy_lammps::error(std::string const &message)
void colvarproxy_lammps::fatal_error(std::string const &message)
{
log(message);
// if (!cvm::debug())
// log("If this error message is unclear, try recompiling the "
// "colvars library and LAMMPS with -DCOLVARS_DEBUG.\n");
_lmp->error->one(FLERR,
"Fatal error in the collective variables module.\n");
}

7
src/USER-COLVARS/colvarproxy_lammps.h
View File

@ -7,10 +7,11 @@
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#ifndef COLVARPROXY_LAMMPS_H
#define COLVARPROXY_LAMMPS_H
#include "colvarproxy_lammps_version.h"
#include "colvarmodule.h"
#include "colvarproxy.h"
#include "colvarvalue.h"
@ -28,10 +29,6 @@
#include <omp.h>
#endif
#ifndef COLVARPROXY_VERSION
#define COLVARPROXY_VERSION "2017-01-09"
#endif
/* struct for packed data communication of coordinates and forces. */
struct commdata {
int tag,type;

10
src/USER-COLVARS/colvarproxy_lammps_version.h Normal file
View File

@ -0,0 +1,10 @@
#ifndef COLVARPROXY_VERSION
#define COLVARPROXY_VERSION "2017-07-15"
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/colvars/colvars
// Please update all Colvars source files before making any changes.
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#endif