Merge pull request #197 from giacomofiorin/colvars_2016-10-05

      for (size_t i = 0; i < colvars.size(); i++) {

        // get total forces (lagging by 1 timestep) from colvars

        // and subtract previous ABF force

        system_force[i] = colvars[i]->total_force().real_value

                        - colvar_forces[i].real_value;

//         if (cvm::debug())

//           cvm::log("ABF System force calc: cv " + cvm::to_str(i) +

//                   " fs " + cvm::to_str(system_force[i]) +

//                   " = ft " + cvm::to_str(colvars[i]->total_force().real_value) +

//                   " - fa " + cvm::to_str(colvar_forces[i].real_value));

        update_system_force(i);

      }

      gradients->acc_force(force_bin, system_force);

    }

        for (size_t i = 0; i < colvars.size(); i++) {

          // If we are outside the range of xi, the force has not been obtained above

          // the function is just an accessor, so cheap to call again anyway

          system_force[i] = colvars[i]->total_force().real_value

                          - colvar_forces[i].real_value;

          update_system_force(i);

        }

        z_gradients->acc_force(z_bin, system_force);

      }

  /// n-dim grid contining CZAR estimator of "real" free energy gradients

  colvar_grid_gradient  *czar_gradients;

  inline int update_system_force(size_t i)

  {

    if (colvars[i]->is_enabled(f_cv_subtract_applied_force)) {

      // this colvar is already subtracting the ABF force

      system_force[i] = colvars[i]->total_force().real_value;

    } else {

      system_force[i] = colvars[i]->total_force().real_value

        - colvar_forces[i].real_value;

    }

    if (cvm::debug())

      cvm::log("ABF System force calc: cv " + cvm::to_str(i) +

               " fs " + cvm::to_str(system_force[i]) +

               " = ft " + cvm::to_str(colvars[i]->total_force().real_value) +

               " - fa " + cvm::to_str(colvar_forces[i].real_value));

    return COLVARS_OK;

  }

  // shared ABF

  bool     shared_on;

  size_t   shared_freq;

std::istream & colvarmodule::read_restart(std::istream &is)

{

  bool warn_total_forces = false;

  {

    // read global restart information

    std::string restart_conf;

      if (restart_version.size() && (restart_version != std::string(COLVARS_VERSION))) {

        cvm::log("This state file was generated with version "+restart_version+"\n");

      }

      if ((restart_version.size() == 0) || (restart_version.compare(std::string(COLVARS_VERSION)) < 0)) {

        // check for total force change

        if (proxy->total_forces_enabled()) {

          warn_total_forces = true;

        }

      }

    }

    is.clear();

    parse->clear_keyword_registry();

  }

  cvm::decrease_depth();

  if (warn_total_forces) {

    cvm::log(cvm::line_marker);

    cvm::log("WARNING:\n");

    std::string const warning("### CHANGES IN THE DEFINITION OF SYSTEM FORCES (NOW TOTAL FORCES)\n\

\n\

Starting from the version 2016-08-10 of the Colvars module, \n\

the role of system forces has been replaced by total forces.\n\

\n\

These include *all* forces acting on a collective variable, whether they\n\

come from the force field potential or from external terms\n\

(e.g. restraints), including forces applied by Colvars itself.\n\

\n\

In NAMD, forces applied by Colvars, IMD, SMD, TMD, symmetry\n\

restraints and tclForces are now all counted in the total force.\n\

\n\

In LAMMPS, forces applied by Colvars itself are now counted in the total\n\

force (all forces from other fixes were being counted already).\n\

\n\

\n\

### WHEN IS THIS CHANGE RELEVANT\n\

\n\

This change affects results *only* when (1) outputSystemForce is\n\

requested or (2) the ABF bias is used.  All other usage cases are\n\

*unaffected* (colvar restraints, metadynamics, etc).\n\

\n\

When system forces are reported (flag: outputSystemForce), their values\n\

in the output may change, but the physical trajectory is never affected.\n\

The physical results of ABF calculations may be affected in some cases.\n\

\n\

\n\

### CHANGES TO ABF CALCULATIONS\n\

\n\

Compared to previous Colvars versions, the ABF method will now attempt\n\

to cancel external forces (for example, boundary walls) and it may be\n\

not possible to resume through a state file a simulation that was\n\

performed with a previous version.\n\

\n\

There are three possible scenarios:\n\

\n\

1. No external forces are applied to the atoms used by ABF: results are\n\

unchanged.\n\

\n\

2. Some of the atoms used by ABF experience external forces, but these\n\

forces are not applied directly to the variables used by ABF\n\

(e.g. another colvar that uses the same atoms, tclForces, etc): in this\n\

case, we recommend beginning a new simulation.\n\

\n\

3. External forces are applied to one or more of the colvars used by\n\

ABF, but no other forces are applied to their atoms: you may use the\n\

subtractAppliedForce keyword inside the corresponding colvars to\n\

continue the previous simulation.\n\n");

    cvm::log(warning);

    cvm::log(cvm::line_marker);

    // update this ahead of time in this special case

    output_prefix = proxy->output_prefix();

    cvm::log("All output files will now be saved with the prefix \""+output_prefix+".tmp.*\".\n");

    output_prefix = output_prefix+".tmp";

    write_output_files();

    cvm::log(cvm::line_marker);

    cvm::log("Please review the important warning above. After that, you may rename:\n\

\""+output_prefix+".tmp.colvars.state\"\n\

to:\n\

\""+output_prefix+".colvars.state\"\n");

    cvm::error("Exiting with error until issue is addressed.\n", FATAL_ERROR);

  }

  return is;

}

#define COLVARMODULE_H

#ifndef COLVARS_VERSION

#define COLVARS_VERSION "2016-09-30"

#define COLVARS_VERSION "2016-10-05"

#endif

#ifndef COLVARS_DEBUG

                 COLVARS_NOT_IMPLEMENTED);

  }

  /// Are total forces being used?

  virtual bool total_forces_enabled() const

  {

    cvm::error("Error: total forces are currently not implemented.\n",

               COLVARS_NOT_IMPLEMENTED);

    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;