Warn that extraBonds angles are cosine-based and allow harmonic 63/4763/3
authorJim Phillips <jim@ks.uiuc.edu>
Thu, 1 Nov 2018 03:40:58 +0000 (22:40 -0500)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Mon, 5 Nov 2018 23:23:07 +0000 (17:23 -0600)
The "normal" flag for extraBonds angles was left uninitialized
when the cosAngles feature was implemented, which tends to
result in extraBonds angles being treated as cosine-based.

Since workflows may unknowingly depend on the current behavior
warn the user and add extraBondsCosAngles option to change it.
Also fix bad comment that reversed meaning of normal flag.

Change-Id: I1a300d6b12654c676c91f15236b135ef3e5fc4e3

src/Molecule.C
src/Parameters.h
src/SimParameters.C
src/SimParameters.h
ug/ug_forcefield.tex

index c5d7e49..c245acd 100644 (file)
@@ -8337,6 +8337,7 @@ void Molecule::build_extra_bonds(Parameters *parameters, StringList *file) {
 //to load
   char err_msg[512];
   int a1,a2,a3,a4; float k, ref;
+  int anglesNormal = ( simParams->extraBondsCosAngles ? 0 : 1 );
   #ifndef MEM_OPT_VERSION
   ResizeArray<Bond> bonds;
   ResizeArray<Angle> angles;
@@ -8413,6 +8414,7 @@ void Molecule::build_extra_bonds(Parameters *parameters, StringList *file) {
         AngleValue tmpv;
         tmpv.k = k;  tmpv.theta0 = ref / 180. * PI;
         tmpv.k_ub = 0;  tmpv.r_ub = 0;
+        tmpv.normal = anglesNormal;
         angle_params.add(tmpv);      
               
       } else if ( ! strncasecmp(type,"dihedral",4) ) {
@@ -8515,6 +8517,14 @@ void Molecule::build_extra_bonds(Parameters *parameters, StringList *file) {
   int extraNumAngles = angle_params.size();
   if ( extraNumAngles ) {
     iout << iINFO << "READ " << extraNumAngles << " EXTRA ANGLES\n" << endi;
+    if ( anglesNormal ) {
+      iout << iINFO << "EXTRA ANGLES ARE NORMAL HARMONIC\n" << endi;
+    } else if ( simParams->extraBondsCosAnglesSetByUser ) {
+      iout << iINFO << "EXTRA ANGLES ARE COSINE-BASED\n" << endi;
+    } else {
+      iout << iWARN << "EXTRA ANGLES ARE COSINE-BASED BY DEFAULT TO MATCH PREVIOUS VERSIONS\n";
+      iout << iWARN << "FOR NORMAL HARMONIC EXTRA ANGLES SET extraBondsCosAngles off\n" << endi;
+    }
     #ifndef MEM_OPT_VERSION
     Angle *newangles = new Angle[numAngles+extraNumAngles];
     memcpy(newangles, this->angles, numAngles*sizeof(Angle));
index c58bb07..3f7b9b7 100644 (file)
@@ -90,7 +90,7 @@ public:
        Real theta0;    //  Rest angle for angle
        Real k_ub;      //  Urey-Bradley force constant
        Real r_ub;      //  Urey-Bradley distance
-  int normal; // Whether we use harmonic (0) or cos-based (1) angle terms
+  int normal; // Whether we use harmonic (1) or cos-based (0) angle terms
 };
 
 typedef struct four_body_consts
index e281623..671ff89 100644 (file)
@@ -2148,6 +2148,9 @@ void SimParameters::config_parser_boundary(ParseOptions &opts) {
    opts.optional("extraBonds", "extraBondsFile",
                "file with list of extra bonds",
                 PARSE_MULTIPLES);
+   opts.optionalB("extraBonds", "extraBondsCosAngles",
+               "Should extra angles be cosine-based to match ancient bug",
+               &extraBondsCosAngles, TRUE);
 
 }
 
@@ -4094,6 +4097,12 @@ void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&
    if ( gridforceOn || mgridforceOn ) {
      parse_mgrid_params(config);
    }
+
+   if ( extraBondsOn ) {
+     extraBondsCosAnglesSetByUser = ! ! config->find("extraBondsCosAngles");
+   } else {
+     extraBondsCosAnglesSetByUser = false;
+   }
       
    if (!opts.defined("constraints"))
    {
index c160da1..4954da7 100644 (file)
@@ -866,6 +866,8 @@ public:
         zVector stirPivot;             // Pivot point of stir axis
 
        Bool extraBondsOn;              // read extra bonded forces
+       Bool extraBondsCosAngles;       // extra angles are cosine-based
+       Bool extraBondsCosAnglesSetByUser; // did the user set this explicitly
 
        Bool consForceOn;               //  Should constant force be applied
   char consForceFile[128];
index 735281a..e5dd11e 100644 (file)
@@ -1254,6 +1254,14 @@ $U(x) = k (1 + cos(n x - x_{ref}))$.
 The only difference between dihedrals and
 impropers is the output field that their potential energy is added to.
 
+Due to a very old bug all NAMD releases prior to 2.13 have used the
+MARTINI cosine-based angle potential function for all extra angles.
+Since workflows may unknowingly depend on this undocumented behavior,
+cosine-based angles remain the default, but a warning is printed
+unless the desired behavior is specified via the new option
+extraBondsCosAngles (defaults to ``on'', set to ``off'' to use
+the normal harmonic angle potential function for all extra angles).
+
 The extra bonded term implementation shares the parallel implementation
 of regular bonded terms in \NAMD, allowing large numbers of extra terms
 to be specified with minimal impact on parallel scalability.
@@ -1270,6 +1278,11 @@ Extra bonded terms are enabled via the following options:
 \NAMDCONFWDEF{extraBonds}{enable extra bonded terms?}{{\tt on} or {\tt off}}{{\tt off}}
 {Specifies whether or not extra bonded terms are present.} 
 
+\item
+\NAMDCONFWDEF{extraBondsCosAngles}{are extra angles cosine-based?}{{\tt on} or {\tt off}}{{\tt on}}
+{Specifies whether or not all extra angles are cosine-based for consistency with previous versions.
+Set to {\tt off} to use normal harmonic angle potential for all extra angles.} 
+
 \item
 \NAMDCONF{extraBondsFile}{file containing extra bonded terms}{file}
 {File containing extra bonded terms.  May be repeated for multiple files.}