Skip to content
Snippets Groups Projects
Commit a3b62b93 authored by Gareth Tribello's avatar Gareth Tribello
Browse files

Added new multicolvar for Angles

parent 647fcee7
No related branches found
No related tags found
No related merge requests found
include ../scripts/test.make
#! FIELDS time a1.mean a2.mean a2.less_than a2.between-1 a2.between-2 a2.more_than
0.000000 1.66115 1.66608 3.40427 0.54041 3.28252 2.45980
0.050000 1.65922 1.66308 3.27471 0.53982 3.15893 2.36446
0.100000 1.65828 1.66185 3.31055 0.53871 3.21053 2.40988
0.150000 1.65822 1.66250 3.41233 0.54348 3.32747 2.47687
0.200000 1.65684 1.66061 3.54666 0.54361 3.45293 2.56758
type=driver
# this is to test a different name
arg="--plumed plumed.dat --trajectory-stride 10 --timestep 0.005 --ixyz trajectory.xyz --dump-forces forces --dump-forces-fmt=%8.4f"
extra_files="../trajectories/trajectory.xyz"
#! FIELDS time parameter a1.less_than a1num.less_than
0.000000 0 0.20260 0.20260
0.000000 1 0.09542 0.09542
0.000000 2 0.09570 0.09570
0.000000 3 -0.09852 -0.10511
0.000000 4 -0.09089 0.01171
0.000000 5 -0.00888 -0.10146
0.000000 6 -0.10407 -0.09748
0.000000 7 -0.00453 -0.10713
0.000000 8 -0.08682 0.00577
0.000000 9 0.18403 0.18403
0.000000 10 0.08177 0.08177
0.000000 11 0.09108 0.09108
0.000000 12 0.08177 0.08177
0.000000 13 0.09142 0.09142
0.000000 14 -0.00615 -0.00615
0.000000 15 0.09108 0.09108
0.000000 16 -0.00615 -0.00615
0.000000 17 0.08454 0.08454
0.050000 0 0.15386 0.15386
0.050000 1 0.07072 0.07072
0.050000 2 0.07001 0.07001
0.050000 3 -0.07411 -0.08090
0.050000 4 -0.06696 0.00765
0.050000 5 -0.00892 -0.07060
0.050000 6 -0.07976 -0.07296
0.050000 7 -0.00376 -0.07837
0.050000 8 -0.06109 0.00059
0.050000 9 0.14840 0.14840
0.050000 10 0.06246 0.06246
0.050000 11 0.07185 0.07185
0.050000 12 0.06246 0.06246
0.050000 13 0.06783 0.06783
0.050000 14 -0.00107 -0.00107
0.050000 15 0.07185 0.07185
0.050000 16 -0.00107 -0.00107
0.050000 17 0.05839 0.05839
0.100000 0 0.11473 0.11473
0.100000 1 0.05034 0.05034
0.100000 2 0.04989 0.04989
0.100000 3 -0.05590 -0.05966
0.100000 4 -0.04763 0.00510
0.100000 5 -0.00691 -0.04948
0.100000 6 -0.05884 -0.05507
0.100000 7 -0.00271 -0.05544
0.100000 8 -0.04297 -0.00041
0.100000 9 0.11709 0.11709
0.100000 10 0.04597 0.04597
0.100000 11 0.05498 0.05498
0.100000 12 0.04597 0.04597
0.100000 13 0.04685 0.04685
0.100000 14 -0.00008 -0.00008
0.100000 15 0.05498 0.05498
0.100000 16 -0.00008 -0.00008
0.100000 17 0.04213 0.04213
0.150000 0 0.08007 0.08007
0.150000 1 0.03377 0.03377
0.150000 2 0.03214 0.03214
0.150000 3 -0.03999 -0.04053
0.150000 4 -0.03153 0.00263
0.150000 5 -0.00383 -0.03271
0.150000 6 -0.04008 -0.03954
0.150000 7 -0.00224 -0.03640
0.150000 8 -0.02831 0.00057
0.150000 9 0.08664 0.08664
0.150000 10 0.03261 0.03261
0.150000 11 0.03806 0.03806
0.150000 12 0.03261 0.03261
0.150000 13 0.02982 0.02982
0.150000 14 -0.00031 -0.00031
0.150000 15 0.03806 0.03806
0.150000 16 -0.00031 -0.00031
0.150000 17 0.02891 0.02891
0.200000 0 0.07023 0.07023
0.200000 1 0.02806 0.02806
0.200000 2 0.02496 0.02496
0.200000 3 -0.03526 -0.03532
0.200000 4 -0.02626 0.00220
0.200000 5 -0.00164 -0.02715
0.200000 6 -0.03497 -0.03491
0.200000 7 -0.00180 -0.03026
0.200000 8 -0.02332 0.00219
0.200000 9 0.07870 0.07870
0.200000 10 0.02833 0.02833
0.200000 11 0.03075 0.03075
0.200000 12 0.02833 0.02833
0.200000 13 0.02444 0.02444
0.200000 14 -0.00168 -0.00168
0.200000 15 0.03075 0.03075
0.200000 16 -0.00168 -0.00168
0.200000 17 0.02343 0.02343
This diff is collapsed.
ANGLES ...
GROUPA=1 GROUPB=2-100
SWITCH={RATIONAL R_0=1.0} MEAN LABEL=a1
... ANGLES
ANGLES ...
GROUPA=1 GROUPB=2-100
SWITCH={RATIONAL R_0=1.0}
MEAN
LESS_THAN={GAUSSIAN R_0=0.5pi}
MORE_THAN={GAUSSIAN R_0=0.5pi}
BETWEEN1={GAUSSIAN LOWER=0.25pi UPPER=0.75pi NORM}
BETWEEN2={GAUSSIAN LOWER=0.25pi UPPER=0.75pi}
NL_STRIDE=20 TOL=0.001
LABEL=a2
... ANGLES
PRINT ARG=a1.*,a2.* FILE=colvar FMT=%8.5f
DUMPDERIVATIVES ARG=a2.* FILE=derivatives FMT=%8.5f
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Copyright (c) 2012 The plumed team
(see the PEOPLE file at the root of the distribution for a list of names)
See http://www.plumed-code.org for more information.
This file is part of plumed, version 2.0.
plumed is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
plumed is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with plumed. If not, see <http://www.gnu.org/licenses/>.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
#include "MultiColvar.h"
#include "tools/Angle.h"
#include "tools/SwitchingFunction.h"
#include "core/ActionRegister.h"
#include <string>
#include <cmath>
using namespace std;
namespace PLMD{
namespace multicolvar{
//+PLUMEDOC MCOLVAR ANGLES
/*
Calculate functions of the distribution of angles .
You can use this command to calculate functions such as:
\f[
f(x) = \sum_{ijk} g( \theta_{ijk} )
\f]
Alternatively you can use this command to calculate functions such as:
\f[
f(x) = \sum_{ijk} s(r_{ij})s(r_{jk}) g(\theta_{ijk})
\f]
where \f$s(r)\f$ is a \ref switchingfunction. This second form means that you can
use this to calculate functions of the angles in the first coordination sphere of
an atom / molecule \cite lj-recon.
\par Examples
The following example instructs plumed to find the average of two angles and to
print it to a file
\verbatim
ANGLES ATOMS=1,2,3 ATOMS=4,5,6 MEAN LABEL=a1
PRINT ARG=a1.mean FILE=colvar
\endverbatim
The following example tells plumed to calculate all angles involving
at least one atom from GROUPA and two atoms from GROUPB in which the distances
are less than 1.0. The number of angles between \f$\frac{\pi}{4}\f$ and
\f$\frac{3\pi}{4}\f$ is then output
\verbatim
ANGLES GROUPA=1-10 GROUPB=11-100 BETWEEN={GAUSSIAN LOWER=0.25pi UPPER=0.75pi} SWITCH={GAUSSIAN R_0=1.0} LABEL=a1
PRINT ARG=a1.between FILE=colvar
\endverbatim
This final example instructs plumed to calculate all the angles in the first coordination
spheres of the atoms. A discretized-normalized histogram of the distribution is then output
\verbatim
ANGLE GROUP=1-38 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=pi NBINS=20} SWTICH=GAUSSIAN R_0=1.0} LABEL=a1
PRINT ARG=a1.* FILE=colvar
\endverbatim
*/
//+ENDPLUMEDOC
class Angles : public MultiColvar {
private:
bool use_sf;
Vector dij, dik;
SwitchingFunction sf1;
SwitchingFunction sf2;
public:
static void registerKeywords( Keywords& keys );
Angles(const ActionOptions&);
// active methods:
virtual double compute( const unsigned& j );
/// Returns the number of coordinates of the field
unsigned getNumberOfFieldDerivatives();
bool contributionIsSmall();
bool isPeriodic(){ return false; }
};
PLUMED_REGISTER_ACTION(Angles,"ANGLES")
void Angles::registerKeywords( Keywords& keys ){
MultiColvar::registerKeywords( keys );
keys.use("ATOMS"); keys.use("MEAN"); keys.use("LESS_THAN");
keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MORE_THAN");
// Could also add Region here in theory
keys.add("atoms-1","GROUP","Calculate angles for each distinct set of three atoms in the group");
keys.add("atoms-2","GROUPA","A group of central atoms about which angles should be calculated");
keys.add("atoms-2","GROUPB","When used in conjuction with GROUPA this keyword instructs plumed "
"to calculate all distinct angles involving one atom from GROUPA "
"and two atoms from GROUPB. The atom from GROUPA is the central atom.");
keys.add("atoms-3","GROUPC","This must be used in conjuction with GROUPA and GROUPB. All angles "
"involving one atom from GROUPA, one atom from GROUPB and one atom from "
"GROUPC are calculated. The GROUPA atoms are assumed to be the central "
"atoms");
keys.add("optional","SWITCH","A switching function that ensures that only angles between atoms that "
"are within a certain fixed cutoff are calculated. The following provides "
"information on the \\ref switchingfunction that are available.");
keys.add("optional","SWITCHA","A switching function on the distance between the atoms in group A and the atoms in "
"group B");
keys.add("optional","SWITCHB","A switching function on the distance between the atoms in group A and the atoms in "
"group B");
}
Angles::Angles(const ActionOptions&ao):
PLUMED_MULTICOLVAR_INIT(ao),
use_sf(false)
{
// Read in the atoms
int natoms=3; readAtoms( natoms );
std::string sfinput,errors; parse("SWITCH",sfinput);
if( sfinput.length()>0 ){
use_sf=true;
weightHasDerivatives=true;
sf1.set(sfinput,errors);
sf2.set(sfinput,errors);
log.printf(" only calculating angles for atoms separated by less than %s\n", sf1.description().c_str() );
} else {
parse("SWITCHA",sfinput);
if(sfinput.length()>0){
use_sf=true;
weightHasDerivatives=true;
sf1.set(sfinput,errors);
sfinput.clear(); parse("SWITCHB",sfinput);
if(sfinput.length()==0) error("found SWITCHA keyword without SWITCHB");
sf2.set(sfinput,errors);
log.printf(" only calculating angles when the distance between GROUPA and GROUPB atoms is less than %s\n", sf1.description().c_str() );
log.printf(" only calculating angles when the distance between GROUPA and GROUPC atoms is less than %s\n", sf2.description().c_str() );
}
}
// And setup the ActionWithVessel
readVesselKeywords();
// And check everything has been read in correctly
checkRead();
}
unsigned Angles::getNumberOfFieldDerivatives(){
return 3*getNumberOfAtoms() + 9;
}
bool Angles::contributionIsSmall(){
dij=getSeparation( getPosition(1), getPosition(2) );
dik=getSeparation( getPosition(1), getPosition(0) );
if(!use_sf) return false;
double w1, w2, dw1, dw2, wtot;
w1=sf1.calculate( dij.modulo(), dw1 );
w2=sf2.calculate( dik.modulo(), dw2 );
wtot=w1*w2; dw1*=w2; dw2*=w1;
if( wtot<getTolerance() ) return true;
addAtomsDerivativeOfWeight( 0, dw2*dik );
addAtomsDerivativeOfWeight( 1, -dw1*dij - dw2*dik );
addAtomsDerivativeOfWeight( 2, dw1*dij );
addBoxDerivativesOfWeight( (-dw1)*Tensor(dij,dij) + (-dw2)*Tensor(dik,dik) );
setWeight( wtot );
return false;
}
double Angles::compute( const unsigned& j ){
Vector ddij,ddik; PLMD::Angle a;
double angle=a.compute(dij,dik,ddij,ddik);
// And finish the calculation
addAtomsDerivatives( 0, ddik );
addAtomsDerivatives( 1, - ddik - ddij );
addAtomsDerivatives( 2, ddij );
addBoxDerivatives( -(Tensor(dij,ddij)+Tensor(dik,ddik)) );
return angle;
}
}
}
......@@ -51,6 +51,17 @@ eprint = {http://www.pnas.org/content/109/14/5196.full.pdf+html},
journal = {Proceedings of the National Academy of Sciences}
}
@article{lj-recon,
title = {Exploring the free energy surfaces of clusters using reconnaissance metadynamics},
publisher = {AIP},
year = {2011},
volume = {135},
number = {11},
pages = {114109},
author = {Gareth A. Tribello and Jérôme Cuny and Hagai Eshet and Michele Parrinello},
journal = { J. Chem. Phys. }
}
@article{camilloni_protG,
Author = {Carlo Camilloni and Davide Provasi and Guido Tiana and Ricardo A. Broglia},
Journal = {Proteins},
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment