Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "tools/SwitchingFunction.h"
23 : #include "core/ActionRegister.h"
24 : #include "multicolvar/AtomValuePack.h"
25 : #include "VectorMultiColvar.h"
26 :
27 : //+PLUMEDOC MCOLVAR BOND_DIRECTIONS
28 : /*
29 : Calculate the vectors connecting atoms that are within cutoff defined using a switching function.
30 :
31 : \par Examples
32 :
33 : */
34 : //+ENDPLUMEDOC
35 :
36 : namespace PLMD {
37 : namespace crystallization {
38 :
39 : class BondOrientation : public VectorMultiColvar {
40 : private:
41 : double rcut2;
42 : SwitchingFunction switchingFunction;
43 : public:
44 : static void registerKeywords( Keywords& keys );
45 : explicit BondOrientation( const ActionOptions& ao );
46 : double calculateWeight( const unsigned& current, const double& weight, multicolvar::AtomValuePack& myatoms ) const override;
47 : void calculateVector( multicolvar::AtomValuePack& myatoms ) const override;
48 : };
49 :
50 10423 : PLUMED_REGISTER_ACTION(BondOrientation,"BOND_DIRECTIONS")
51 :
52 3 : void BondOrientation::registerKeywords( Keywords& keys ) {
53 3 : VectorMultiColvar::registerKeywords( keys );
54 6 : keys.add("numbered","ATOMS","the atoms involved in each of the vectors you wish to calculate. "
55 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one vector will be "
56 : "calculated for each ATOM keyword you specify (all ATOM keywords should "
57 : "specify the indices of two atoms). The eventual number of quantities calculated by this "
58 : "action will depend on what functions of the distribution you choose to calculate.");
59 6 : keys.reset_style("ATOMS","atoms");
60 6 : keys.add("atoms-1","GROUP","Calculate the distance between each distinct pair of atoms in the group");
61 6 : keys.add("atoms-2","GROUPA","Calculate the distances between all the atoms in GROUPA and all "
62 : "the atoms in GROUPB. This must be used in conjunction with GROUPB.");
63 6 : keys.add("atoms-2","GROUPB","Calculate the distances between all the atoms in GROUPA and all the atoms "
64 : "in GROUPB. This must be used in conjunction with GROUPA.");
65 6 : keys.add("compulsory","NN","12","The n parameter of the switching function ");
66 6 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
67 6 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
68 6 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
69 6 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
70 : "The following provides information on the \\ref switchingfunction that are available. "
71 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
72 6 : keys.use("VMEAN"); keys.use("VSUM");
73 3 : }
74 :
75 2 : BondOrientation::BondOrientation( const ActionOptions& ao ):
76 : Action(ao),
77 2 : VectorMultiColvar(ao)
78 : {
79 : // Read atoms
80 2 : weightHasDerivatives=true;
81 : std::vector<AtomNumber> all_atoms;
82 4 : readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
83 2 : if( atom_lab.size()==0 ) readAtomsLikeKeyword( "ATOMS", 2, all_atoms );
84 2 : setupMultiColvarBase( all_atoms );
85 : // Read in the switching function
86 4 : std::string sw, errors; parse("SWITCH",sw);
87 2 : if(sw.length()>0) {
88 2 : switchingFunction.set(sw,errors);
89 : } else {
90 0 : double r_0=-1.0, d_0; int nn, mm;
91 0 : parse("NN",nn); parse("MM",mm);
92 0 : parse("R_0",r_0); parse("D_0",d_0);
93 0 : if( r_0<0.0 ) error("you must set a value for R_0");
94 0 : switchingFunction.set(nn,mm,r_0,d_0);
95 : }
96 2 : log.printf(" orientation of those bonds with lengths are less than %s\n",( switchingFunction.description() ).c_str() );
97 : // Set the link cell cutoff
98 2 : setLinkCellCutoff( switchingFunction.get_dmax() );
99 2 : double rcut = switchingFunction.get_dmax(); rcut2 = rcut*rcut;
100 : // Set the dimensionality of the vectors
101 2 : setVectorDimensionality(3);
102 2 : }
103 :
104 150 : double BondOrientation::calculateWeight( const unsigned& current, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
105 150 : Vector distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) ); double distm=distance.modulo2();
106 150 : if( distm>rcut2 ) return 0.0;
107 60 : double df, ww=switchingFunction.calculateSqr( distm, df );
108 : // Derivatives of weights
109 60 : addAtomDerivatives( 0, 0, -df*weight*distance, myatoms );
110 60 : addAtomDerivatives( 0, 1, df*weight*distance, myatoms );
111 60 : myatoms.addBoxDerivatives( 0, (-df)*weight*Tensor(distance,distance) );
112 60 : return ww;
113 : }
114 :
115 60 : void BondOrientation::calculateVector( multicolvar::AtomValuePack& myatoms ) const {
116 60 : Vector distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
117 :
118 60 : addAtomDerivatives( 2, 0, Vector(-1.0,0,0), myatoms );
119 60 : addAtomDerivatives( 2, 1, Vector(+1.0,0,0), myatoms );
120 60 : myatoms.addBoxDerivatives( 2, Tensor(distance,Vector(-1.0,0,0)) );
121 60 : myatoms.addValue( 2, distance[0] );
122 :
123 60 : addAtomDerivatives( 3, 0, Vector(0,-1.0,0), myatoms );
124 60 : addAtomDerivatives( 3, 1, Vector(0,+1.0,0), myatoms );
125 60 : myatoms.addBoxDerivatives( 3, Tensor(distance,Vector(0,-1.0,0)) );
126 60 : myatoms.addValue( 3, distance[1] );
127 :
128 60 : addAtomDerivatives( 4, 0, Vector(0,0,-1.0), myatoms );
129 60 : addAtomDerivatives( 4, 1, Vector(0,0,+1.0), myatoms );
130 60 : myatoms.addBoxDerivatives( 4, Tensor(distance,Vector(0,0,-1.0)) );
131 60 : myatoms.addValue( 4, distance[2] );
132 60 : }
133 :
134 : }
135 : }
|
