diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf index 5a6b253363..8fcc37c26e 100644 Binary files a/doc/src/PDF/colvars-refman-lammps.pdf and b/doc/src/PDF/colvars-refman-lammps.pdf differ diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index 8d54ce755b..159e39d075 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -41,7 +41,7 @@ Syntax * template-ID(post-reacted) = ID of a molecule template containing post-reaction topology * map_file = name of file specifying corresponding atom-IDs in the pre- and post-reacted templates * zero or more individual keyword/value pairs may be appended to each react argument -* individual_keyword = *prob* or *max_rxn* or *stabilize_steps* or *update_edges* +* individual_keyword = *prob* or *max_rxn* or *stabilize_steps* or *custom_charges* .. parsed-literal:: @@ -52,10 +52,9 @@ Syntax N = maximum number of reactions allowed to occur *stabilize_steps* value = timesteps timesteps = number of timesteps to apply the internally-created :doc:`nve/limit ` fix to reacting atoms - *update_edges* value = *none* or *charges* or *custom* - *none* = do not update topology near the edges of reaction templates - *charges* = update atomic charges of all atoms in reaction templates - *custom* = force the update of user-specified atomic charges + *custom_charges* value = *no* or *fragmentID* + no = update all atomic charges (default) + fragmentID = ID of molecule fragment whose charges are updated Examples """""""" @@ -271,14 +270,14 @@ A discussion of correctly handling this is also provided on the The map file is a text document with the following format: A map file has a header and a body. The header of map file the -contains one mandatory keyword and five optional keywords. The +contains one mandatory keyword and four optional keywords. The mandatory keyword is 'equivalences': .. parsed-literal:: N *equivalences* = # of atoms N in the reaction molecule templates -The optional keywords are 'edgeIDs', 'deleteIDs', 'customIDs' and +The optional keywords are 'edgeIDs', 'deleteIDs', 'chiralIDs' and 'constraints': .. parsed-literal:: @@ -286,10 +285,9 @@ The optional keywords are 'edgeIDs', 'deleteIDs', 'customIDs' and N *edgeIDs* = # of edge atoms N in the pre-reacted molecule template N *deleteIDs* = # of atoms N that are specified for deletion N *chiralIDs* = # of specified chiral centers N - N *customIDs* = # of atoms N that are specified for a custom update N *constraints* = # of specified reaction constraints N -The body of the map file contains two mandatory sections and five +The body of the map file contains two mandatory sections and four optional sections. The first mandatory section begins with the keyword 'BondingIDs' and lists the atom IDs of the bonding atom pair in the pre-reacted molecule template. The second mandatory section begins @@ -303,16 +301,11 @@ molecule template. The second optional section begins with the keyword 'DeleteIDs' and lists the atom IDs of pre-reaction template atoms to delete. The third optional section begins with the keyword 'ChiralIDs' lists the atom IDs of chiral atoms whose handedness should be -enforced. The fourth optional section begins with the keyword 'Custom -Edges' and allows for forcing the update of a specific atom's atomic -charge. The first column is the ID of an atom near the edge of the -pre-reacted molecule template, and the value of the second column is -either 'none' or 'charges.' Further details are provided in the -discussion of the 'update_edges' keyword. The fifth optional section -begins with the keyword 'Constraints' and lists additional criteria -that must be satisfied in order for the reaction to occur. Currently, -there are five types of constraints available, as discussed below: -'distance', 'angle', 'dihedral', 'arrhenius', and 'rmsd'. +enforced. The fourth optional section begins with the keyword +'Constraints' and lists additional criteria that must be satisfied in +order for the reaction to occur. Currently, there are five types of +constraints available, as discussed below: 'distance', 'angle', +'dihedral', 'arrhenius', and 'rmsd'. A sample map file is given below: @@ -488,17 +481,12 @@ individually tuned for each fix reaction step. Note that in some situations, decreasing rather than increasing this parameter will result in an increase in stability. -The *update_edges* keyword can increase the number of atoms whose -atomic charges are updated, when the pre-reaction template contains -edge atoms. When the value is set to 'charges,' all atoms' atomic -charges are updated to those specified by the post-reaction template, -including atoms near the edge of reaction templates. When the value is -set to 'custom,' an additional section must be included in the map -file that specifies whether or not to update charges, on a per-atom -basis. The format of this section is detailed above. Listing a -pre-reaction atom ID with a value of 'charges' will force the update -of the atom's charge, even if it is near a template edge. Atoms not -near a template edge are unaffected by this setting. +The *custom_charges* keyword can be used to specify which atoms' +atomic charges are updated. When the value is set to 'no,' all atomic +charges are updated to those specified by the post-reaction template +(default). Otherwise, the value should be the name of a molecule +fragment defined in the pre-reaction molecule template. In this case, +only the atomic charges of atoms in the molecule fragment are updated. A few other considerations: @@ -584,7 +572,7 @@ Default """"""" The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60, -reset_mol_ids = yes, update_edges = none +reset_mol_ids = yes, custom_charges = no ---------- diff --git a/lib/colvars/Makefile.deps b/lib/colvars/Makefile.deps index ef7cb9bef0..78e8768edd 100644 --- a/lib/colvars/Makefile.deps +++ b/lib/colvars/Makefile.deps @@ -4,79 +4,227 @@ $(COLVARS_OBJ_DIR)colvaratoms.o: colvaratoms.cpp colvarmodule.h \ colvarproxy_tcl.h colvarproxy_volmaps.h colvarparse.h colvarparams.h \ colvaratoms.h colvardeps.h $(COLVARS_OBJ_DIR)colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h \ - colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h \ - colvarproxy_tcl.h colvarproxy_volmaps.h colvar.h colvarparse.h \ - colvarparams.h colvardeps.h colvarbias_abf.h colvarbias.h colvargrid.h \ - colvar_UIestimator.h + colvars_version.h colvar.h colvarvalue.h colvartypes.h colvarparse.h \ + colvarparams.h colvardeps.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarbias_abf.h colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ + colvarbias.h colvargrid.h colvar_UIestimator.h $(COLVARS_OBJ_DIR)colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h \ colvars_version.h colvarbias.h colvar.h colvarvalue.h colvartypes.h \ - colvarparse.h colvarparams.h colvardeps.h colvarbias_alb.h + colvarparse.h colvarparams.h colvardeps.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarbias_alb.h $(COLVARS_OBJ_DIR)colvarbias.o: colvarbias.cpp colvarmodule.h \ colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h \ colvarproxy_tcl.h colvarproxy_volmaps.h colvarbias.h colvar.h \ - colvarparse.h colvarparams.h colvardeps.h colvargrid.h + colvarparse.h colvarparams.h colvardeps.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvargrid.h $(COLVARS_OBJ_DIR)colvarbias_histogram.o: colvarbias_histogram.cpp \ colvarmodule.h colvars_version.h colvarproxy.h colvartypes.h \ colvarvalue.h colvarproxy_tcl.h colvarproxy_volmaps.h colvar.h \ - colvarparse.h colvarparams.h colvardeps.h colvarbias_histogram.h \ - colvarbias.h colvargrid.h + colvarparse.h colvarparams.h colvardeps.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarbias_histogram.h colvarbias.h colvargrid.h $(COLVARS_OBJ_DIR)colvarbias_meta.o: colvarbias_meta.cpp colvarmodule.h \ colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h \ colvarproxy_tcl.h colvarproxy_volmaps.h colvar.h colvarparse.h \ - colvarparams.h colvardeps.h colvarbias_meta.h colvarbias.h colvargrid.h + colvarparams.h colvardeps.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarbias_meta.h colvarbias.h colvargrid.h $(COLVARS_OBJ_DIR)colvarbias_restraint.o: colvarbias_restraint.cpp \ colvarmodule.h colvars_version.h colvarproxy.h colvartypes.h \ colvarvalue.h colvarproxy_tcl.h colvarproxy_volmaps.h \ colvarbias_restraint.h colvarbias.h colvar.h colvarparse.h \ - colvarparams.h colvardeps.h + colvarparams.h colvardeps.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h $(COLVARS_OBJ_DIR)colvarcomp_angles.o: colvarcomp_angles.cpp \ colvarmodule.h colvars_version.h colvar.h colvarvalue.h colvartypes.h \ - colvarparse.h colvarparams.h colvardeps.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparse.h colvarparams.h colvardeps.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp_apath.o: colvarcomp_apath.cpp colvarmodule.h \ colvars_version.h colvarvalue.h colvartypes.h colvarparse.h \ - colvarparams.h colvar.h colvardeps.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparams.h colvar.h colvardeps.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp_coordnums.o: colvarcomp_coordnums.cpp \ colvarmodule.h colvars_version.h colvarparse.h colvarvalue.h \ colvartypes.h colvarparams.h colvaratoms.h colvarproxy.h \ colvarproxy_tcl.h colvarproxy_volmaps.h colvardeps.h colvar.h \ + lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ colvarcomp.h colvar_arithmeticpath.h colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp.o: colvarcomp.cpp colvarmodule.h \ colvars_version.h colvarvalue.h colvartypes.h colvar.h colvarparse.h \ - colvarparams.h colvardeps.h colvarcomp.h colvaratoms.h colvarproxy.h \ - colvarproxy_tcl.h colvarproxy_volmaps.h colvar_arithmeticpath.h \ - colvar_geometricpath.h + colvarparams.h colvardeps.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp_distances.o: colvarcomp_distances.cpp \ colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \ - colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \ - colvaratoms.h colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparse.h colvarparams.h colvar.h colvardeps.h \ + lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp_gpath.o: colvarcomp_gpath.cpp colvarmodule.h \ colvars_version.h colvarvalue.h colvartypes.h colvarparse.h \ - colvarparams.h colvar.h colvardeps.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparams.h colvar.h colvardeps.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp_protein.o: colvarcomp_protein.cpp \ colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \ - colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \ - colvaratoms.h colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparse.h colvarparams.h colvar.h colvardeps.h \ + lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp_rotations.o: colvarcomp_rotations.cpp \ colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \ - colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \ - colvaratoms.h colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparse.h colvarparams.h colvar.h colvardeps.h \ + lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp_volmaps.o: colvarcomp_volmaps.cpp \ colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \ - colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \ - colvaratoms.h colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparse.h colvarparams.h colvar.h colvardeps.h \ + lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvar.o: colvar.cpp colvarmodule.h colvars_version.h \ colvarvalue.h colvartypes.h colvarparse.h colvarparams.h colvar.h \ - colvardeps.h colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_tcl.h \ + colvardeps.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_tcl.h \ colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h \ colvarscript.h colvarbias.h colvarscript_commands.h \ colvarscript_commands_colvar.h colvarscript_commands_bias.h @@ -86,18 +234,38 @@ $(COLVARS_OBJ_DIR)colvardeps.o: colvardeps.cpp colvarmodule.h \ colvarparams.h $(COLVARS_OBJ_DIR)colvargrid.o: colvargrid.cpp colvarmodule.h \ colvars_version.h colvarvalue.h colvartypes.h colvarparse.h \ - colvarparams.h colvar.h colvardeps.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h colvargrid.h + colvarparams.h colvar.h colvardeps.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h \ + colvargrid.h $(COLVARS_OBJ_DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h \ colvars_version.h colvarparse.h colvarvalue.h colvartypes.h \ colvarparams.h colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar.h colvardeps.h colvarbias.h colvarbias_abf.h colvargrid.h \ - colvar_UIestimator.h colvarbias_alb.h colvarbias_histogram.h \ - colvarbias_meta.h colvarbias_restraint.h colvarscript.h \ - colvarscript_commands.h colvarscript_commands_colvar.h \ - colvarscript_commands_bias.h colvaratoms.h colvarcomp.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvar.h colvardeps.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarbias.h colvarbias_abf.h colvargrid.h colvar_UIestimator.h \ + colvarbias_alb.h colvarbias_histogram.h colvarbias_meta.h \ + colvarbias_restraint.h colvarscript.h colvarscript_commands.h \ + colvarscript_commands_colvar.h colvarscript_commands_bias.h \ + colvaratoms.h colvarcomp.h colvar_arithmeticpath.h \ + colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarparams.o: colvarparams.cpp colvarmodule.h \ colvars_version.h colvarvalue.h colvartypes.h colvarparams.h $(COLVARS_OBJ_DIR)colvarparse.o: colvarparse.cpp colvarmodule.h \ @@ -107,6 +275,15 @@ $(COLVARS_OBJ_DIR)colvarproxy.o: colvarproxy.cpp colvarmodule.h \ colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h \ colvarproxy_tcl.h colvarproxy_volmaps.h colvarscript.h colvarbias.h \ colvar.h colvarparse.h colvarparams.h colvardeps.h \ + lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ colvarscript_commands.h colvarscript_commands_colvar.h \ colvarscript_commands_bias.h colvaratoms.h $(COLVARS_OBJ_DIR)colvarproxy_replicas.o: colvarproxy_replicas.cpp \ @@ -122,28 +299,66 @@ $(COLVARS_OBJ_DIR)colvarscript.o: colvarscript.cpp colvarproxy.h \ colvarmodule.h colvars_version.h colvartypes.h colvarvalue.h \ colvarproxy_tcl.h colvarproxy_volmaps.h colvardeps.h colvarparse.h \ colvarparams.h colvarscript.h colvarbias.h colvar.h \ + lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ colvarscript_commands.h colvarscript_commands_colvar.h \ colvarscript_commands_bias.h $(COLVARS_OBJ_DIR)colvarscript_commands.o: colvarscript_commands.cpp \ colvarproxy.h colvarmodule.h colvars_version.h colvartypes.h \ colvarvalue.h colvarproxy_tcl.h colvarproxy_volmaps.h colvardeps.h \ colvarparse.h colvarparams.h colvarscript.h colvarbias.h colvar.h \ + lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ colvarscript_commands.h colvarscript_commands_colvar.h \ colvarscript_commands_bias.h $(COLVARS_OBJ_DIR)colvarscript_commands_bias.o: \ colvarscript_commands_bias.cpp colvarproxy.h colvarmodule.h \ colvars_version.h colvartypes.h colvarvalue.h colvarproxy_tcl.h \ colvarproxy_volmaps.h colvardeps.h colvarparse.h colvarparams.h \ - colvarscript.h colvarbias.h colvar.h colvarscript_commands.h \ - colvarscript_commands_colvar.h colvarscript_commands_bias.h + colvarscript.h colvarbias.h colvar.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarscript_commands.h colvarscript_commands_colvar.h \ + colvarscript_commands_bias.h $(COLVARS_OBJ_DIR)colvarscript_commands_colvar.o: \ colvarscript_commands_colvar.cpp colvarproxy.h colvarmodule.h \ colvars_version.h colvartypes.h colvarvalue.h colvarproxy_tcl.h \ colvarproxy_volmaps.h colvardeps.h colvarparse.h colvarparams.h \ - colvarscript.h colvarbias.h colvar.h colvarscript_commands.h \ - colvarscript_commands_colvar.h colvarscript_commands_bias.h + colvarscript.h colvarbias.h colvar.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarscript_commands.h colvarscript_commands_colvar.h \ + colvarscript_commands_bias.h $(COLVARS_OBJ_DIR)colvartypes.o: colvartypes.cpp colvarmodule.h \ colvars_version.h colvartypes.h colvarparse.h colvarvalue.h \ - colvarparams.h + colvarparams.h ../../src/math_eigen.h $(COLVARS_OBJ_DIR)colvarvalue.o: colvarvalue.cpp colvarmodule.h \ colvars_version.h colvarvalue.h colvartypes.h diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index fc8b490631..ed6fec1ac5 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -18,7 +18,9 @@ #include "colvarcomp.h" #include "colvarscript.h" - +#if (__cplusplus >= 201103L) +std::map> colvar::global_cvc_map = std::map>(); +#endif colvar::colvar() { @@ -724,15 +726,15 @@ int colvar::init_output_flags(std::string const &conf) return COLVARS_OK; } - - - // read the configuration and set up corresponding instances, for // each type of component implemented template int colvar::init_components_type(std::string const &conf, char const * /* def_desc */, char const *def_config_key) { +#if (__cplusplus >= 201103L) + global_cvc_map[def_config_key] = [](const std::string& cvc_conf){return new def_class_name(cvc_conf);}; +#endif size_t def_count = 0; std::string def_conf = ""; size_t pos = 0; diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index 98905dbe17..5cf0077a39 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -12,6 +12,11 @@ #include +#if (__cplusplus >= 201103L) +#include +#include +#endif + #include "colvarmodule.h" #include "colvarvalue.h" #include "colvarparse.h" @@ -114,6 +119,7 @@ public: /// List of biases that depend on this colvar std::vector biases; + protected: @@ -605,6 +611,14 @@ public: // components that do not handle any atoms directly class map_total; + /// getter of the global cvc map +#if (__cplusplus >= 201103L) + /// A global mapping of cvc names to the cvc constructors + static const std::map>& get_global_cvc_map() { + return global_cvc_map; + } +#endif + protected: /// \brief Array of \link colvar::cvc \endlink objects @@ -639,6 +653,11 @@ protected: double dev_null; #endif +#if (__cplusplus >= 201103L) + /// A global mapping of cvc names to the cvc constructors + static std::map> global_cvc_map; +#endif + public: /// \brief Sorted array of (zero-based) IDs for all atoms involved diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index 785e4c6edc..57df68ab6e 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -1493,8 +1493,6 @@ class colvar::linearCombination : public colvar::cvc { protected: - /// Map from string to the types of colvar components - std::map> string_cv_map; /// Sub-colvar components std::vector cv; /// If all sub-cvs use explicit gradients then we also use it @@ -1514,8 +1512,6 @@ class colvar::CVBasedPath : public colvar::cvc { protected: - /// Map from string to the types of colvar components - std::map> string_cv_map; /// Sub-colvar components std::vector cv; /// Reference colvar values from path diff --git a/lib/colvars/colvarcomp_gpath.cpp b/lib/colvars/colvarcomp_gpath.cpp index d563403628..8849f9dbe4 100644 --- a/lib/colvars/colvarcomp_gpath.cpp +++ b/lib/colvars/colvarcomp_gpath.cpp @@ -19,10 +19,6 @@ #include "colvar.h" #include "colvarcomp.h" -namespace GeometricPathCV { -void init_string_cv_map(std::map>& string_cv_map); -} - bool compareColvarComponent(colvar::cvc *i, colvar::cvc *j) { return i->name < j->name; @@ -406,9 +402,8 @@ void colvar::gzpath::apply_force(colvarvalue const &force) { } colvar::linearCombination::linearCombination(std::string const &conf): cvc(conf) { - GeometricPathCV::init_string_cv_map(string_cv_map); // Lookup all available sub-cvcs - for (auto it_cv_map = string_cv_map.begin(); it_cv_map != string_cv_map.end(); ++it_cv_map) { + for (auto it_cv_map = colvar::get_global_cvc_map().begin(); it_cv_map != colvar::get_global_cvc_map().end(); ++it_cv_map) { if (key_lookup(conf, it_cv_map->first.c_str())) { std::vector sub_cvc_confs; get_key_string_multi_value(conf, it_cv_map->first.c_str(), sub_cvc_confs); @@ -506,9 +501,8 @@ void colvar::linearCombination::apply_force(colvarvalue const &force) { } colvar::CVBasedPath::CVBasedPath(std::string const &conf): cvc(conf) { - GeometricPathCV::init_string_cv_map(string_cv_map); // Lookup all available sub-cvcs - for (auto it_cv_map = string_cv_map.begin(); it_cv_map != string_cv_map.end(); ++it_cv_map) { + for (auto it_cv_map = colvar::get_global_cvc_map().begin(); it_cv_map != colvar::get_global_cvc_map().end(); ++it_cv_map) { if (key_lookup(conf, it_cv_map->first.c_str())) { std::vector sub_cvc_confs; get_key_string_multi_value(conf, it_cv_map->first.c_str(), sub_cvc_confs); @@ -909,33 +903,4 @@ void colvar::gzpathCV::apply_force(colvarvalue const &force) { } } -void GeometricPathCV::init_string_cv_map(std::map>& string_cv_map) { - string_cv_map["distance"] = [](const std::string& conf){return new colvar::distance(conf);}; - string_cv_map["dihedral"] = [](const std::string& conf){return new colvar::dihedral(conf);}; - string_cv_map["angle"] = [](const std::string& conf){return new colvar::angle(conf);}; - string_cv_map["rmsd"] = [](const std::string& conf){return new colvar::rmsd(conf);}; - string_cv_map["gyration"] = [](const std::string& conf){return new colvar::gyration(conf);}; - string_cv_map["inertia"] = [](const std::string& conf){return new colvar::inertia(conf);}; - string_cv_map["inertiaZ"] = [](const std::string& conf){return new colvar::inertia_z(conf);}; - string_cv_map["tilt"] = [](const std::string& conf){return new colvar::tilt(conf);}; - string_cv_map["distanceZ"] = [](const std::string& conf){return new colvar::distance_z(conf);}; - string_cv_map["distanceXY"] = [](const std::string& conf){return new colvar::distance_xy(conf);}; - string_cv_map["polarTheta"] = [](const std::string& conf){return new colvar::polar_theta(conf);}; - string_cv_map["polarPhi"] = [](const std::string& conf){return new colvar::polar_phi(conf);}; - string_cv_map["distanceVec"] = [](const std::string& conf){return new colvar::distance_vec(conf);}; - string_cv_map["orientationAngle"] = [](const std::string& conf){return new colvar::orientation_angle(conf);}; - string_cv_map["distancePairs"] = [](const std::string& conf){return new colvar::distance_pairs(conf);}; - string_cv_map["dipoleMagnitude"] = [](const std::string& conf){return new colvar::dipole_magnitude(conf);}; - string_cv_map["coordNum"] = [](const std::string& conf){return new colvar::coordnum(conf);}; - string_cv_map["selfCoordNum"] = [](const std::string& conf){return new colvar::selfcoordnum(conf);}; - string_cv_map["dipoleAngle"] = [](const std::string& conf){return new colvar::dipole_angle(conf);}; - string_cv_map["orientation"] = [](const std::string& conf){return new colvar::orientation(conf);}; - string_cv_map["orientationProj"] = [](const std::string& conf){return new colvar::orientation_proj(conf);}; - string_cv_map["eigenvector"] = [](const std::string& conf){return new colvar::eigenvector(conf);}; - string_cv_map["cartesian"] = [](const std::string& conf){return new colvar::cartesian(conf);}; - string_cv_map["alpha"] = [](const std::string& conf){return new colvar::alpha_angles(conf);}; - string_cv_map["dihedralPC"] = [](const std::string& conf){return new colvar::dihedPC(conf);}; - string_cv_map["linearCombination"] = [](const std::string& conf){return new colvar::linearCombination(conf);}; -} - #endif diff --git a/lib/colvars/colvars_version.h b/lib/colvars/colvars_version.h index f69d55a818..33e05d72d1 100644 --- a/lib/colvars/colvars_version.h +++ b/lib/colvars/colvars_version.h @@ -1,3 +1,3 @@ #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2020-07-07" +#define COLVARS_VERSION "2020-09-17" #endif diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp index 75d7384a00..3cbaed63a8 100644 --- a/lib/colvars/colvartypes.cpp +++ b/lib/colvars/colvartypes.cpp @@ -14,25 +14,19 @@ #include "colvartypes.h" #include "colvarparse.h" +#ifdef COLVARS_LAMMPS +// Use open-source Jacobi implementation +#include "math_eigen_impl.h" +#else +// Fall back to NR routine +#include "nr_jacobi.h" +#endif + bool colvarmodule::rotation::monitor_crossings = false; cvm::real colvarmodule::rotation::crossing_threshold = 1.0E-02; -namespace { - -/// Numerical recipes diagonalization -static int jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot); - -/// Eigenvector sort -static int eigsrt(cvm::real *d, cvm::real **v); - -/// Transpose the matrix -static int transpose(cvm::real **v); - -} - - std::string cvm::rvector::to_simple_string() const { std::ostringstream os; @@ -248,6 +242,66 @@ cvm::quaternion::position_derivative_inner(cvm::rvector const &pos, // Seok C, Dill KA. Using quaternions to calculate RMSD. J Comput // Chem. 25(15):1849-57 (2004) DOI: 10.1002/jcc.20110 PubMed: 15376254 +#ifdef COLVARS_LAMMPS +namespace { + inline void *new_Jacobi_solver(int size) { + return reinterpret_cast(new MathEigen::Jacobi &, + cvm::matrix2d &>(4)); + } +} +#endif + + +colvarmodule::rotation::rotation() +{ + b_debug_gradients = false; +#ifdef COLVARS_LAMMPS + jacobi = new_Jacobi_solver(4); +#else + jacobi = NULL; +#endif +} + + +colvarmodule::rotation::rotation(cvm::quaternion const &qi) + : q(qi) +{ + b_debug_gradients = false; +#ifdef COLVARS_LAMMPS + jacobi = new_Jacobi_solver(4); +#else + jacobi = NULL; +#endif +} + + +colvarmodule::rotation::rotation(cvm::real angle, cvm::rvector const &axis) +{ + b_debug_gradients = false; + cvm::rvector const axis_n = axis.unit(); + cvm::real const sina = cvm::sin(angle/2.0); + q = cvm::quaternion(cvm::cos(angle/2.0), + sina * axis_n.x, sina * axis_n.y, sina * axis_n.z); +#ifdef COLVARS_LAMMPS + jacobi = new_Jacobi_solver(4); +#else + jacobi = NULL; +#endif +} + + +colvarmodule::rotation::~rotation() +{ +#ifdef COLVARS_LAMMPS + delete reinterpret_cast< + MathEigen::Jacobi &, + cvm::matrix2d &> *>(jacobi); +#endif +} + + void colvarmodule::rotation::build_correlation_matrix( std::vector const &pos1, std::vector const &pos2) @@ -291,10 +345,12 @@ void colvarmodule::rotation::compute_overlap_matrix() } -void colvarmodule::rotation::diagonalize_matrix( - cvm::matrix2d &m, - cvm::vector1d &eigval, - cvm::matrix2d &eigvec) +#ifndef COLVARS_LAMMPS +namespace { + +void diagonalize_matrix(cvm::matrix2d &m, + cvm::vector1d &eigval, + cvm::matrix2d &eigvec) { eigval.resize(4); eigval.reset(); @@ -303,15 +359,15 @@ void colvarmodule::rotation::diagonalize_matrix( // diagonalize int jac_nrot = 0; - if (jacobi(m.c_array(), eigval.c_array(), eigvec.c_array(), &jac_nrot) != + if (NR_Jacobi::jacobi(m.c_array(), eigval.c_array(), eigvec.c_array(), &jac_nrot) != COLVARS_OK) { - cvm::error("Too many iterations in routine jacobi.\n" + cvm::error("Too many iterations in jacobi diagonalization.\n" "This is usually the result of an ill-defined set of atoms for " "rotational alignment (RMSD, rotateReference, etc).\n"); } - eigsrt(eigval.c_array(), eigvec.c_array()); + NR_Jacobi::eigsrt(eigval.c_array(), eigvec.c_array()); // jacobi saves eigenvectors by columns - transpose(eigvec.c_array()); + NR_Jacobi::transpose(eigvec.c_array()); // normalize eigenvectors for (size_t ie = 0; ie < 4; ie++) { @@ -327,6 +383,9 @@ void colvarmodule::rotation::diagonalize_matrix( } } +} +#endif + // Calculate the rotation, plus its derivatives @@ -349,7 +408,28 @@ void colvarmodule::rotation::calc_optimal_rotation( cvm::log("S = "+cvm::to_str(S_backup, cvm::cv_width, cvm::cv_prec)+"\n"); } + S_eigval.resize(4); + S_eigvec.resize(4, 4); + +#ifdef COLVARS_LAMMPS + MathEigen::Jacobi &, + cvm::matrix2d &> *ecalc = + reinterpret_cast &, + cvm::matrix2d &> *>(jacobi); + + int ierror = ecalc->Diagonalize(S, S_eigval, S_eigvec); + if (ierror) { + cvm::error("Too many iterations in jacobi diagonalization.\n" + "This is usually the result of an ill-defined set of atoms for " + "rotational alignment (RMSD, rotateReference, etc).\n"); + } +#else diagonalize_matrix(S, S_eigval, S_eigvec); +#endif + + // eigenvalues and eigenvectors cvm::real const L0 = S_eigval[0]; cvm::real const L1 = S_eigval[1]; @@ -521,7 +601,11 @@ void colvarmodule::rotation::calc_optimal_rotation( // cvm::log("S_new = "+cvm::to_str(cvm::to_str (S_new), cvm::cv_width, cvm::cv_prec)+"\n"); +#ifdef COLVARS_LAMMPS + ecalc->Diagonalize(S_new, S_new_eigval, S_new_eigvec); +#else diagonalize_matrix(S_new, S_new_eigval, S_new_eigvec); +#endif cvm::real const &L0_new = S_new_eigval[0]; cvm::quaternion const Q0_new(S_new_eigvec[0]); @@ -544,138 +628,3 @@ void colvarmodule::rotation::calc_optimal_rotation( -// Numerical Recipes routine for diagonalization - -#define ROTATE(a,i,j,k,l) g=a[i][j]; \ - h=a[k][l]; \ - a[i][j]=g-s*(h+g*tau); \ - a[k][l]=h+s*(g-h*tau); - -#define n 4 - - -namespace { - -int jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot) -{ - int j,iq,ip,i; - cvm::real tresh,theta,tau,t,sm,s,h,g,c; - - cvm::vector1d b(n); - cvm::vector1d z(n); - - for (ip=0;ip 4 && (cvm::real)(cvm::fabs(d[ip])+g) == (cvm::real)cvm::fabs(d[ip]) - && (cvm::real)(cvm::fabs(d[iq])+g) == (cvm::real)cvm::fabs(d[iq])) - a[ip][iq]=0.0; - else if (cvm::fabs(a[ip][iq]) > tresh) { - h=d[iq]-d[ip]; - if ((cvm::real)(cvm::fabs(h)+g) == (cvm::real)cvm::fabs(h)) - t=(a[ip][iq])/h; - else { - theta=0.5*h/(a[ip][iq]); - t=1.0/(cvm::fabs(theta)+cvm::sqrt(1.0+theta*theta)); - if (theta < 0.0) t = -t; - } - c=1.0/cvm::sqrt(1+t*t); - s=t*c; - tau=s/(1.0+c); - h=t*a[ip][iq]; - z[ip] -= h; - z[iq] += h; - d[ip] -= h; - d[iq] += h; - a[ip][iq]=0.0; - for (j=0;j<=ip-1;j++) { - ROTATE(a,j,ip,j,iq) - } - for (j=ip+1;j<=iq-1;j++) { - ROTATE(a,ip,j,j,iq) - } - for (j=iq+1;j= p) p=d[k=j]; - if (k != i) { - d[k]=d[i]; - d[i]=p; - for (j=0;j const &pos2); /// Default constructor - inline rotation() - : b_debug_gradients(false) - {} + rotation(); /// Constructor after a quaternion - inline rotation(cvm::quaternion const &qi) - : q(qi), - b_debug_gradients(false) - { - } + rotation(cvm::quaternion const &qi); /// Constructor after an axis of rotation and an angle (in radians) - inline rotation(cvm::real angle, cvm::rvector const &axis) - : b_debug_gradients(false) - { - cvm::rvector const axis_n = axis.unit(); - cvm::real const sina = cvm::sin(angle/2.0); - q = cvm::quaternion(cvm::cos(angle/2.0), - sina * axis_n.x, sina * axis_n.y, sina * axis_n.z); - } + rotation(cvm::real angle, cvm::rvector const &axis); /// Destructor - inline ~rotation() - {} + ~rotation(); /// Return the rotated vector inline cvm::rvector rotate(cvm::rvector const &v) const @@ -1432,7 +1418,6 @@ public: return q.rotation_matrix(); } - /// \brief Return the spin angle (in degrees) with respect to the /// provided axis (which MUST be normalized) inline cvm::real spin_angle(cvm::rvector const &axis) const @@ -1537,10 +1522,8 @@ protected: /// Compute the overlap matrix S (used by calc_optimal_rotation()) void compute_overlap_matrix(); - /// Diagonalize a given matrix m (used by calc_optimal_rotation()) - static void diagonalize_matrix(cvm::matrix2d &m, - cvm::vector1d &eigval, - cvm::matrix2d &eigvec); + /// Pointer to instance of Jacobi solver + void *jacobi; }; diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp index 9f9a1ce00d..814f8a0fcf 100644 --- a/src/USER-REACTION/fix_bond_react.cpp +++ b/src/USER-REACTION/fix_bond_react.cpp @@ -201,7 +201,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : memory->create(seed,nreacts,"bond/react:seed"); memory->create(limit_duration,nreacts,"bond/react:limit_duration"); memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag"); - memory->create(update_edges_flag,nreacts,"bond/react:update_edges_flag"); + memory->create(custom_charges_fragid,nreacts,"bond/react:custom_charges_fragid"); memory->create(constraints,1,MAXCONARGS,"bond/react:constraints"); memory->create(var_flag,NUMVARVALS,nreacts,"bond/react:var_flag"); memory->create(var_id,NUMVARVALS,nreacts,"bond/react:var_id"); @@ -221,7 +221,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : seed[i] = 12345; max_rxn[i] = INT_MAX; stabilize_steps_flag[i] = 0; - update_edges_flag[i] = 0; + custom_charges_fragid[i] = -1; // set default limit duration to 60 timesteps limit_duration[i] = 60; reaction_count[i] = 0; @@ -367,13 +367,15 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : limit_duration[rxn] = utils::numeric(FLERR,arg[iarg+1],false,lmp); stabilize_steps_flag[rxn] = 1; iarg += 2; - } else if (strcmp(arg[iarg],"update_edges") == 0) { + } else if (strcmp(arg[iarg],"custom_charges") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: " - "'update_edges' has too few arguments"); - if (strcmp(arg[iarg+1],"none") == 0) update_edges_flag[rxn] = 0; - else if (strcmp(arg[iarg+1],"charges") == 0) update_edges_flag[rxn] = 1; - else if (strcmp(arg[iarg+1],"custom") == 0) update_edges_flag[rxn] = 2; - else error->all(FLERR,"Illegal value for 'update_edges' keyword'"); + "'custom_charges' has too few arguments"); + if (strcmp(arg[iarg+1],"no") == 0) custom_charges_fragid[rxn] = -1; //default + else { + custom_charges_fragid[rxn] = atom->molecules[unreacted_mol[rxn]]->findfragment(arg[iarg+1]); + if (custom_charges_fragid[rxn] < 0) error->one(FLERR,"Bond/react: Molecule fragment for " + "'custom_charges' keyword does not exist"); + } iarg += 2; } else error->all(FLERR,"Illegal fix bond/react command: unknown keyword"); } @@ -389,15 +391,14 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : memory->create(reverse_equiv,max_natoms,2,nreacts,"bond/react:reverse_equiv"); memory->create(edge,max_natoms,nreacts,"bond/react:edge"); memory->create(landlocked_atoms,max_natoms,nreacts,"bond/react:landlocked_atoms"); - memory->create(custom_edges,max_natoms,nreacts,"bond/react:custom_edges"); + memory->create(custom_charges,max_natoms,nreacts,"bond/react:custom_charges"); memory->create(delete_atoms,max_natoms,nreacts,"bond/react:delete_atoms"); memory->create(chiral_atoms,max_natoms,6,nreacts,"bond/react:chiral_atoms"); for (int j = 0; j < nreacts; j++) for (int i = 0; i < max_natoms; i++) { edge[i][j] = 0; - if (update_edges_flag[j] == 1) custom_edges[i][j] = 1; - else custom_edges[i][j] = 0; + custom_charges[i][j] = 1; // update all partial charges by default delete_atoms[i][j] = 0; for (int k = 0; k < 6; k++) { chiral_atoms[i][k][j] = 0; @@ -419,6 +420,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : iatomtype[i] = onemol->type[ibonding[i]-1]; jatomtype[i] = onemol->type[jbonding[i]-1]; find_landlocked_atoms(i); + if (custom_charges_fragid[i] >= 0) CustomCharges(custom_charges_fragid[i],i); } // initialize Marsaglia RNG with processor-unique seed (Arrhenius prob) @@ -436,7 +438,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : } delete [] files; - if (atom->molecular != Atom::ATOMIC) + if (atom->molecular != Atom::MOLECULAR) error->all(FLERR,"Bond/react: Cannot use fix bond/react with non-molecular systems"); // check if bonding atoms are 1-2, 1-3, or 1-4 bonded neighbors @@ -528,7 +530,7 @@ FixBondReact::~FixBondReact() memory->destroy(equivalences); memory->destroy(reverse_equiv); memory->destroy(landlocked_atoms); - memory->destroy(custom_edges); + memory->destroy(custom_charges); memory->destroy(delete_atoms); memory->destroy(chiral_atoms); @@ -545,7 +547,7 @@ FixBondReact::~FixBondReact() memory->destroy(var_flag); memory->destroy(var_id); memory->destroy(stabilize_steps_flag); - memory->destroy(update_edges_flag); + memory->destroy(custom_charges_fragid); memory->destroy(iatomtype); memory->destroy(jatomtype); @@ -2641,11 +2643,11 @@ void FixBondReact::update_everything() twomol = atom->molecules[reacted_mol[rxnID]]; for (int j = 0; j < twomol->natoms; j++) { int jj = equivalences[j][1][rxnID]-1; - if ((landlocked_atoms[j][rxnID] == 1 || custom_edges[jj][rxnID] == 1) && - atom->map(update_mega_glove[jj+1][i]) >= 0 && + if (atom->map(update_mega_glove[jj+1][i]) >= 0 && atom->map(update_mega_glove[jj+1][i]) < nlocal) { - type[atom->map(update_mega_glove[jj+1][i])] = twomol->type[j]; - if (twomol->qflag && atom->q_flag) { + if (landlocked_atoms[j][rxnID] == 1) + type[atom->map(update_mega_glove[jj+1][i])] = twomol->type[j]; + if (twomol->qflag && atom->q_flag && custom_charges[jj][rxnID] == 1) { double *q = atom->q; q[atom->map(update_mega_glove[jj+1][i])] = twomol->q[j]; } @@ -3155,7 +3157,6 @@ void FixBondReact::read(int myrxn) error->one(FLERR,"Bond/react: Number of equivalences in map file must " "equal number of atoms in reaction templates"); } - else if (strstr(line,"customIDs")) sscanf(line,"%d",&ncustom); else if (strstr(line,"deleteIDs")) sscanf(line,"%d",&ndelete); else if (strstr(line,"chiralIDs")) sscanf(line,"%d",&nchiral); else if (strstr(line,"constraints")) { @@ -3171,7 +3172,7 @@ void FixBondReact::read(int myrxn) // loop over sections of superimpose file - int equivflag = 0, bondflag = 0, customedgesflag = 0; + int equivflag = 0, bondflag = 0; while (strlen(keyword)) { if (strcmp(keyword,"BondingIDs") == 0) { bondflag = 1; @@ -3188,9 +3189,6 @@ void FixBondReact::read(int myrxn) } else if (strcmp(keyword,"Equivalences") == 0) { equivflag = 1; Equivalences(line, myrxn); - } else if (strcmp(keyword,"Custom Edges") == 0) { - customedgesflag = 1; - CustomEdges(line, myrxn); } else if (strcmp(keyword,"DeleteIDs") == 0) { DeleteAtoms(line, myrxn); } else if (strcmp(keyword,"ChiralIDs") == 0) { @@ -3206,12 +3204,6 @@ void FixBondReact::read(int myrxn) // error check if (bondflag == 0 || equivflag == 0) error->all(FLERR,"Bond/react: Map file missing BondingIDs or Equivalences section\n"); - - if (update_edges_flag[myrxn] == 2 && customedgesflag == 0) - error->all(FLERR,"Bond/react: Map file must have a Custom Edges section when using 'update_edges custom'\n"); - - if (update_edges_flag[myrxn] != 2 && customedgesflag == 1) - error->all(FLERR,"Bond/react: Specify 'update_edges custom' to include Custom Edges section in map file\n"); } void FixBondReact::EdgeIDs(char *line, int myrxn) @@ -3246,28 +3238,6 @@ void FixBondReact::Equivalences(char *line, int myrxn) } } -void FixBondReact::CustomEdges(char *line, int myrxn) -{ - // 0 for 'none', 1 for 'charges' - - int tmp; - int n = MAX(strlen("none"),strlen("charges")) + 1; - char *edgemode = new char[n]; - for (int i = 0; i < ncustom; i++) { - readline(line); - sscanf(line,"%d %s",&tmp,edgemode); - if (tmp > onemol->natoms) - error->one(FLERR,"Bond/react: Invalid template atom ID in map file"); - if (strcmp(edgemode,"none") == 0) - custom_edges[tmp-1][myrxn] = 0; - else if (strcmp(edgemode,"charges") == 0) - custom_edges[tmp-1][myrxn] = 1; - else - error->one(FLERR,"Bond/react: Illegal value in 'Custom Edges' section of map file"); - } - delete [] edgemode; -} - void FixBondReact::DeleteAtoms(char *line, int myrxn) { int tmp; @@ -3280,6 +3250,15 @@ void FixBondReact::DeleteAtoms(char *line, int myrxn) } } +void FixBondReact::CustomCharges(int ifragment, int myrxn) +{ + for (int i = 0; i < onemol->natoms; i++) + if (onemol->fragmentmask[ifragment][i]) + custom_charges[i][myrxn] = 1; + else + custom_charges[i][myrxn] = 0; +} + void FixBondReact::ChiralCenters(char *line, int myrxn) { int tmp; diff --git a/src/USER-REACTION/fix_bond_react.h b/src/USER-REACTION/fix_bond_react.h index 61a1bd4213..895a63b7ae 100644 --- a/src/USER-REACTION/fix_bond_react.h +++ b/src/USER-REACTION/fix_bond_react.h @@ -64,7 +64,7 @@ class FixBondReact : public Fix { int reset_mol_ids_flag; int custom_exclude_flag; int *stabilize_steps_flag; - int *update_edges_flag; + int *custom_charges_fragid; int nconstraints; int narrhenius; double **constraints; @@ -113,7 +113,7 @@ class FixBondReact : public Fix { int *ibonding,*jbonding; int *closeneigh; // indicates if bonding atoms of a rxn are 1-2, 1-3, or 1-4 neighbors - int nedge,nequivalent,ncustom,ndelete,nchiral,nconstr; // # edge, equivalent, custom atoms in mapping file + int nedge,nequivalent,ndelete,nchiral,nconstr; // # edge, equivalent atoms in mapping file int attempted_rxn; // there was an attempt! int *local_rxn_count; int *ghostly_rxn_count; @@ -127,7 +127,7 @@ class FixBondReact : public Fix { int ***equivalences; // relation between pre- and post-reacted templates int ***reverse_equiv; // re-ordered equivalences int **landlocked_atoms; // all atoms at least three bonds away from edge atoms - int **custom_edges; // atoms in molecule templates with incorrect valences + int **custom_charges; // atoms whose charge should be updated int **delete_atoms; // atoms in pre-reacted templates to delete int ***chiral_atoms; // pre-react chiral atoms. 1) flag 2) orientation 3-4) ordered atom types @@ -151,8 +151,8 @@ class FixBondReact : public Fix { void read(int); void EdgeIDs(char *, int); void Equivalences(char *, int); - void CustomEdges(char *, int); void DeleteAtoms(char *, int); + void CustomCharges(int, int); void ChiralCenters(char *, int); void Constraints(char *, int); void readID(char *, int, int, int);