Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 : Copyright (c) 2012-2017 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 "core/ActionShortcut.h" 23 : #include "core/ActionRegister.h" 24 : #include "core/PlumedMain.h" 25 : #include "core/ActionSet.h" 26 : #include "core/ActionWithValue.h" 27 : #include "CoordinationNumbers.h" 28 : #include "multicolvar/MultiColvarShortcuts.h" 29 : 30 : #include <complex> 31 : 32 : namespace PLMD { 33 : namespace symfunc { 34 : 35 : //+PLUMEDOC MCOLVAR LOCAL_CRYSTALINITY 36 : /* 37 : Calculate the local crystalinity symmetry function 38 : 39 : \par Examples 40 : 41 : 42 : */ 43 : //+ENDPLUMEDOC 44 : 45 : 46 : class LocalCrystallinity : public ActionShortcut { 47 : public: 48 : static void registerKeywords( Keywords& keys ); 49 : explicit LocalCrystallinity(const ActionOptions&); 50 : }; 51 : 52 : PLUMED_REGISTER_ACTION(LocalCrystallinity,"LOCAL_CRYSTALINITY") 53 : 54 3 : void LocalCrystallinity::registerKeywords( Keywords& keys ) { 55 3 : CoordinationNumbers::shortcutKeywords( keys ); 56 6 : keys.add("numbered","GVECTOR","the coefficients of the linear combinations to compute for the CV"); 57 6 : keys.needsAction("ONES"); keys.needsAction("MATRIX_VECTOR_PRODUCT"); 58 6 : keys.needsAction("COMBINE"); keys.needsAction("CUSTOM"); 59 3 : } 60 : 61 1 : LocalCrystallinity::LocalCrystallinity( const ActionOptions& ao): 62 : Action(ao), 63 1 : ActionShortcut(ao) 64 : { 65 : // This builds an adjacency matrix 66 3 : std::string sp_str, specA, specB; parse("SPECIES",sp_str); parse("SPECIESA",specA); parse("SPECIESB",specB); 67 1 : CoordinationNumbers::expandMatrix( true, getShortcutLabel(), sp_str, specA, specB, this ); 68 : // Input for denominator (coord) 69 1 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_mat"); 70 1 : plumed_assert( av && av->getNumberOfComponents()>0 && (av->copyOutput(0))->getRank()==2 ); 71 1 : std::string size; Tools::convert( (av->copyOutput(0))->getShape()[1], size ); 72 2 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + size ); 73 2 : readInputLine( getShortcutLabel() + "_denom: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_mat.w," + getShortcutLabel() + "_ones"); 74 : // Input for numerator 75 1 : std::string finput = ""; 76 1 : for(unsigned i=1;; ++i) { 77 4 : std::vector<std::string> gvec; std::string istr; Tools::convert( i, istr ); 78 8 : if( !parseNumberedVector("GVECTOR",i,gvec) ) { break; } 79 3 : if( gvec.size()!=3 ) error("gvectors should have size 3"); 80 : // This is the dot product between the input gvector and the bond vector 81 9 : readInputLine( getShortcutLabel() + "_dot" + istr + ": COMBINE ARG=" + getShortcutLabel() + "_mat.x," + getShortcutLabel() + "_mat.y," + getShortcutLabel() + "_mat.z " 82 6 : "PERIODIC=NO COEFFICIENTS=" + gvec[0] + "," + gvec[1] + "," + gvec[2] ); 83 : // Now calculate the sine and cosine of the dot product 84 6 : readInputLine( getShortcutLabel() + "_cos" + istr + ": CUSTOM ARG=" + getShortcutLabel() +"_mat.w," + getShortcutLabel() + "_dot" + istr + " FUNC=x*cos(y) PERIODIC=NO"); 85 6 : readInputLine( getShortcutLabel() + "_sin" + istr + ": CUSTOM ARG=" + getShortcutLabel() +"_mat.w," + getShortcutLabel() + "_dot" + istr + " FUNC=x*sin(y) PERIODIC=NO"); 86 : // And sum up the sine and cosine over the coordination spheres 87 6 : readInputLine( getShortcutLabel() + "_cossum" + istr + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_cos" + istr + "," + getShortcutLabel() + "_ones"); 88 6 : readInputLine( getShortcutLabel() + "_sinsum" + istr + ": MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_sin" + istr + "," + getShortcutLabel() + "_ones"); 89 : // And average the sine and cosine over the number of bonds 90 6 : readInputLine( getShortcutLabel() + "_cosmean" + istr + ": CUSTOM FUNC=x/y PERIODIC=NO ARG=" + getShortcutLabel() + "_cossum" + istr + "," + getShortcutLabel() + "_denom"); 91 6 : readInputLine( getShortcutLabel() + "_sinmean" + istr + ": CUSTOM FUNC=x/y PERIODIC=NO ARG=" + getShortcutLabel() + "_sinsum" + istr + "," + getShortcutLabel() + "_denom"); 92 : // And work out the square modulus of this complex number 93 6 : readInputLine( getShortcutLabel() + "_" + istr + ": CUSTOM FUNC=x*x+y*y PERIODIC=NO ARG=" + getShortcutLabel() + "_cosmean" + istr + "," + getShortcutLabel() + "_sinmean" + istr); 94 : // These are all the kvectors that we are adding together in the final combine for the final CV 95 3 : if( i>1 ) finput += ","; 96 6 : finput += getShortcutLabel() + "_" + istr; 97 4 : } 98 : // This computes the final CV 99 2 : readInputLine( getShortcutLabel() + ": COMBINE NORMALIZE PERIODIC=NO ARG=" + finput ); 100 : // Now calculate the total length of the vector 101 2 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel(), getShortcutLabel(), "", this ); 102 1 : } 103 : 104 : } 105 : } 106 :