Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2014-2020 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             : #ifndef __PLUMED_symfunc_Fccubic_h
      23             : #define __PLUMED_symfunc_Fccubic_h
      24             : #include "function/FunctionSetup.h"
      25             : 
      26             : #include <string>
      27             : #include <cmath>
      28             : 
      29             : namespace PLMD {
      30             : namespace symfunc {
      31             : 
      32             : template<typename T=double>
      33             : struct Fccubic {
      34             :   using precision=T;
      35             :   precision alpha;
      36             :   precision a1;
      37             :   precision b1;
      38             :   static void registerKeywords( Keywords& keys );
      39             :   static void read( Fccubic& func,
      40             :                     ActionWithArguments* action,
      41             :                     function::FunctionOptions& funcout );
      42             :   static void calc( const Fccubic& func,
      43             :                     bool noderiv,
      44             :                     View<const precision> args,
      45             :                     function::FunctionOutput& funcout );
      46             : 
      47             : #ifdef __PLUMED_HAS_OPENACC
      48             :   void toACCDevice() const {
      49             : #pragma acc enter data copyin(this[0:1], alpha,a1,b1)
      50             :   }
      51             :   void removeFromACCDevice() const {
      52             : #pragma acc exit data delete(b1,a1,alpha,this[0:1])
      53             :   }
      54             : #endif // __PLUMED_HAS_OPENACC
      55             : };
      56             : 
      57             : template <typename T>
      58          32 : void Fccubic<T>::registerKeywords( Keywords& keys ) {
      59          32 :   keys.add("compulsory","ALPHA","3.0",
      60             :            "The alpha parameter of the angular function");
      61          64 :   keys.setValueDescription("matrix",
      62             :                            "a function that measures the similarity with an fcc environment");
      63          32 : }
      64             : 
      65             : template <typename T>
      66           6 : void Fccubic<T>::read( Fccubic& func,
      67             :                        ActionWithArguments* action,
      68             :                        function::FunctionOptions& funcout ) {
      69             :   // Scaling factors such that '1' corresponds to fcc lattice
      70             :   // and '0' corresponds to isotropic (liquid)
      71           6 :   action->parse("ALPHA",func.alpha);
      72           6 :   func.a1 = 80080. / (2717. + 16*func.alpha);
      73           6 :   func.b1 = 16.*(func.alpha-143)/(2717+16*func.alpha);
      74           6 :   action->log.printf("  setting alpha paramter equal to %f \n",func.alpha);
      75           6 : }
      76             : 
      77             : template <typename T>
      78      328574 : void Fccubic<T>::calc( const Fccubic& func,
      79             :                        bool noderiv,
      80             :                        const View<const precision> args,
      81             :                        function::FunctionOutput& funcout ) {
      82      328574 :   precision x2 = args[0]*args[0];
      83      328574 :   precision x4 = x2*x2;
      84             : 
      85      328574 :   precision y2 = args[1]*args[1];
      86      328574 :   precision y4 = y2*y2;
      87             : 
      88      328574 :   precision z2 = args[2]*args[2];
      89      328574 :   precision z4 = z2*z2;
      90             : 
      91      328574 :   precision d2 = x2 + y2 + z2;
      92      328574 :   if( d2 < epsilon ) {
      93             :     d2 = 1;
      94             :   }
      95      328574 :   precision r8 = pow( d2, 4 );
      96      328574 :   precision r12 = pow( d2, 6 );
      97             : 
      98      328574 :   precision tmp = ((x4*y4)+(x4*z4)+(y4*z4))/r8-func.alpha*x4*y4*z4/r12;
      99             : 
     100      328574 :   precision t0 = (x2*y4+x2*z4)/r8-func.alpha*x2*y4*z4/r12;
     101      328574 :   precision t1 = (y2*x4+y2*z4)/r8-func.alpha*y2*x4*z4/r12;
     102      328574 :   precision t2 = (z2*x4+z2*y4)/r8-func.alpha*z2*x4*y4/r12;
     103      328574 :   precision t3 = (2*tmp-func.alpha*x4*y4*z4/r12)/d2;
     104             : 
     105      328574 :   if( !noderiv ) {
     106        4166 :     funcout.derivs[0][0]=4*func.a1*args[0]*(t0-t3);
     107        4166 :     funcout.derivs[0][1]=4*func.a1*args[1]*(t1-t3);
     108        4166 :     funcout.derivs[0][2]=4*func.a1*args[2]*(t2-t3);
     109             :   }
     110             : 
     111             :   // Set the value and the derivatives
     112      328574 :   funcout.values[0] = (func.a1*tmp+func.b1);
     113      328574 : }
     114             : 
     115             : }
     116             : }
     117             : 
     118             : #endif //__PLUMED_symfunc_Fccubic_h