Logo Search packages:      
Sourcecode: adun.app version File versions  Download package

AdEnzymixNonbonded.c

/*
   Project: Adun

   Copyright (C) 2005 Michael Johnston & Jordi Villa-Freixa

   Author: Michael Johnston

   This application is free software; you can redistribute it and/or
   modify it under the terms of the GNU General Public
   License as published by the Free Software Foundation; either
   version 2 of the License, or (at your option) any later version.

   This application 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
   Library General Public License for more details.

   You should have received a copy of the GNU General Public
   License along with this library; if not, write to the Free
   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111 USA.
*/

#include <Base/AdForceFieldFunctions.h>


inline void AdEnzymixNonbondedForce(ListElement* interaction, double* vdw_pot, double* est_pot, 
                        double EPSILON_RP, double** coordinates, double** accelerations,
                        double cut)
{
      int atom_one, atom_two;
      double force_mag;
      double length_rec, vdw_hold, est_hold;
      double A, B;
      double charge_A, charge_B;
      Vector3D seperation_s;
      
      atom_one = interaction->bond[0];
      atom_two = interaction->bond[1];
                  
      //calculate seperation vectir (r1 - r2)
      
      *(seperation_s.vector + 0) = coordinates[atom_one][0] - coordinates[atom_two][0];
      *(seperation_s.vector + 1) = coordinates[atom_one][1] - coordinates[atom_two][1];
      *(seperation_s.vector + 2) = coordinates[atom_one][2] - coordinates[atom_two][2];
      
      //translate the information in the bond
      
      charge_A = coordinates[atom_one][5];
      charge_B = coordinates[atom_two][5];

      //calculate the interatomic distance

      Ad3DVectorLength(&seperation_s);
      
      //add this length to the linked list wlement the bond belongs to
      
      interaction->length = seperation_s.length;

/*    if(seperation_s.length > cut)
            return;*/
                  
      //get reciprocal of seperation

      length_rec = 1/seperation_s.length;
      
      //calculate vdw holder
      
      vdw_hold = pow(length_rec, 6);

      //while this is being calculate 
      //calculate some non dependant variables

      //est holder

      est_hold = EPSILON_RP*charge_A*charge_B*length_rec;

      //vdw holder may be finished now so calculate vdw consts
      
      A = interaction->params[0]*vdw_hold*vdw_hold;
      B = interaction->params[1]*vdw_hold;
      
      *est_pot += est_hold;
      *vdw_pot += A - B;      

      force_mag = est_hold*length_rec;
      
      //add the vdw force to the est force
      
      force_mag +=  6*length_rec*(2*A - B);
      force_mag *= length_rec;

      //calculate the force on atom one along the vector (r1 - r2)
      //the force on atom two is the opposite of this force

      *(seperation_s.vector + 0) *= force_mag;
      *(seperation_s.vector + 1) *= force_mag;
      *(seperation_s.vector + 2) *= force_mag;
      
      accelerations[atom_one][0] += *(seperation_s.vector + 0);
      accelerations[atom_one][1] += *(seperation_s.vector + 1);
      accelerations[atom_one][2] += *(seperation_s.vector + 2);

      accelerations[atom_two][0] -= *(seperation_s.vector + 0);
      accelerations[atom_two][1] -= *(seperation_s.vector + 1);
      accelerations[atom_two][2] -= *(seperation_s.vector + 2);

}

inline void AdEnzymixShiftedNonbondedForce(ListElement* interaction, double* vdw_pot, double* est_pot, 
                        double EPSILON_RP, double** coordinates, double** accelerations,
                        double cut, float r_cutoff2)
{
      int atom_one, atom_two;
      double force_mag;
      double length_rec, vdw_hold, est_hold;
      double A, B;
      double charge_A, charge_B, shift_fac;
      Vector3D seperation_s;
      
      atom_one = interaction->bond[0];
      atom_two = interaction->bond[1];

      //printf("Cutoff is %lf. Reciprocal of cutoff squared is %f\n", cut,r_cutoff2);
                  
      //calculate seperation vectir (r1 - r2)
      
      *(seperation_s.vector + 0) = coordinates[atom_one][0] - coordinates[atom_two][0];
      *(seperation_s.vector + 1) = coordinates[atom_one][1] - coordinates[atom_two][1];
      *(seperation_s.vector + 2) = coordinates[atom_one][2] - coordinates[atom_two][2];
      
      //translate the information in the bond
      
      charge_A = coordinates[atom_one][5];
      charge_B = coordinates[atom_two][5];

      //calculate the interatomic distance

      Ad3DVectorLength(&seperation_s);
      
      //add this length to the linked list wlement the bond belongs to
      
      interaction->length = seperation_s.length;

/*    if(seperation_s.length > cut)
            return;*/
                  
      //get reciprocal of seperation

      length_rec = 1/seperation_s.length;
      
      //calculate vdw holder
      
      vdw_hold = pow(length_rec, 6);

      //while this is being calculate 
      //calculate some non dependant variables

      //est holder

      est_hold = EPSILON_RP*charge_A*charge_B*length_rec;
      
      //est shift factor
      
      shift_fac = (cut - seperation_s.length)*(cut - seperation_s.length)*r_cutoff2;

      //vdw holder may be finished now so calculate vdw consts
      
      A = interaction->params[0]*vdw_hold*vdw_hold;
      B = interaction->params[1]*vdw_hold;

      //the shifted est force is EST_HOLD*(length_rec - r_cutoff2*length)
      
      force_mag = (length_rec - r_cutoff2*seperation_s.length);

      //caclulate the potentials    

      *est_pot += (est_hold*shift_fac);
      *vdw_pot += A - B;      
      
      //apply the shift to the est force

      force_mag *= est_hold;

      //add the vdw force

      force_mag +=  6*length_rec*(2*A - B);
      force_mag *= length_rec;

      //calculate the force on atom one along the vector (r1 - r2)
      //the force on atom two is the opposite of this force

      *(seperation_s.vector + 0) *= force_mag;
      *(seperation_s.vector + 1) *= force_mag;
      *(seperation_s.vector + 2) *= force_mag;
      
      accelerations[atom_one][0] += *(seperation_s.vector + 0);
      accelerations[atom_one][1] += *(seperation_s.vector + 1);
      accelerations[atom_one][2] += *(seperation_s.vector + 2);

      accelerations[atom_two][0] -= *(seperation_s.vector + 0);
      accelerations[atom_two][1] -= *(seperation_s.vector + 1);
      accelerations[atom_two][2] -= *(seperation_s.vector + 2);
}

inline void AdEnzymixGRFNonbondedForce(ListElement* interaction, double* vdw_pot, double* est_pot, 
                        double EPSILON_RP, double** coordinates, double** accelerations,
                        double cut, double b0, double b1)
{
      int atom_one, atom_two;
      double force_mag;
      double length_rec, vdw_hold, est_hold;
      double A, B;
      double charge_A, charge_B;
      Vector3D seperation_s;
      
      atom_one = interaction->bond[0];
      atom_two = interaction->bond[1];
                  
      //calculate seperation vectir (r1 - r2)
      
      *(seperation_s.vector + 0) = coordinates[atom_one][0] - coordinates[atom_two][0];
      *(seperation_s.vector + 1) = coordinates[atom_one][1] - coordinates[atom_two][1];
      *(seperation_s.vector + 2) = coordinates[atom_one][2] - coordinates[atom_two][2];
      
      //translate the information in the bond
      
      charge_A = coordinates[atom_one][5];
      charge_B = coordinates[atom_two][5];

      //calculate the interatomic distance

      Ad3DVectorLength(&seperation_s);
      
      //add this length to the linked list wlement the bond belongs to
      
      interaction->length = seperation_s.length;

/*    if(seperation_s.length > cut)
            return;*/
                  
      //get reciprocal of seperation

      length_rec = 1/seperation_s.length;
      
      //calculate vdw holder
      
      vdw_hold = pow(length_rec, 6);

      //while this is being calculate 
      //calculate some non dependant variables

      //est holder

      est_hold = EPSILON_RP*charge_A*charge_B*length_rec;

      //vdw holder may be finished now so calculate vdw consts
      
      A = interaction->params[0]*vdw_hold*vdw_hold;
      B = interaction->params[1]*vdw_hold;
      
      *est_pot += est_hold + EPSILON_RP*charge_B*charge_A*b0;
      *vdw_pot += A - B;      
      
      force_mag = est_hold*length_rec + EPSILON_RP*charge_B*charge_A*b1*seperation_s.length;

      //add the vdw force to the est force
      
      force_mag +=  6*length_rec*(2*A - B);
      force_mag *=length_rec;

      //calculate the force on atom one along the vector (r1 - r2)
      //the force on atom two is the opposite of this force

      *(seperation_s.vector + 0) *= force_mag;
      *(seperation_s.vector + 1) *= force_mag;
      *(seperation_s.vector + 2) *= force_mag;
      
      accelerations[atom_one][0] += *(seperation_s.vector + 0);
      accelerations[atom_one][1] += *(seperation_s.vector + 1);
      accelerations[atom_one][2] += *(seperation_s.vector + 2);

      accelerations[atom_two][0] -= *(seperation_s.vector + 0);
      accelerations[atom_two][1] -= *(seperation_s.vector + 1);
      accelerations[atom_two][2] -= *(seperation_s.vector + 2);
}

inline void AdEnzymixNonbondedEnergy(ListElement *Interaction, double *vdw_pot, double *est_pot, double EPSILON_RP, double **coordinates,
                        double cut)
{
      int atom_one, atom_two;
      double length_rec, vdw_hold, est_hold;
      double A, B;
      double charge_A, charge_B;
      Vector3D seperation_s;
      
      atom_one = Interaction->bond[0];
      atom_two = Interaction->bond[1];
                  
      //calculate seperation vectir (r1 - r2)
      
      *(seperation_s.vector + 0) = coordinates[atom_one][0] - coordinates[atom_two][0];
      *(seperation_s.vector + 1) = coordinates[atom_one][1] - coordinates[atom_two][1];
      *(seperation_s.vector + 2) = coordinates[atom_one][2] - coordinates[atom_two][2];
      
      //translate the information in the bond
      
      charge_A = coordinates[atom_one][5];
      charge_B = coordinates[atom_two][5];

      //calculate the interatomic distance

      Ad3DVectorLength(&seperation_s);
      
      //add this length to the linked list wlement the bond belongs to
      
      Interaction->length = seperation_s.length;

/*    if(seperation_s.length > cut)
            return;*/
                  
      //get reciprocal of seperation

      length_rec = 1/seperation_s.length;
      
      //calculate vdw holder
      
      vdw_hold = pow(length_rec, 6);

      //while this is being calculate 
      //calculate some non dependant variables

      //est holder

      est_hold = (EPSILON_RP*charge_A*charge_B*length_rec);

      //vdw holder may be finished now so calculate vdw consts
      
      A = Interaction->params[0]*vdw_hold*vdw_hold;
      B = Interaction->params[1]*vdw_hold;

/*    if(atom_two==0)
            printf("%-12lf%-12lf%-12lf%-12lf\n", charge_A, charge_B, seperation_s.length, est_hold);*/
      
      *est_pot += est_hold;
      *vdw_pot += A - B;      
}

inline void AdEnzymixShiftedNonbondedEnergy(ListElement *Interaction, double *vdw_pot, double *est_pot,  
                  double EPSILON_RP, double **coordinates, double cut, float r_cutoff2)
{           
      int atom_one, atom_two;
      double length_rec, vdw_hold, est_hold;
      double A, B;
      double charge_A, charge_B, shift_fac;
      Vector3D seperation_s;
      
      atom_one = Interaction->bond[0];
      atom_two = Interaction->bond[1];
                  
      //calculate seperation vectir (r1 - r2)
      
      *(seperation_s.vector + 0) = coordinates[atom_one][0] - coordinates[atom_two][0];
      *(seperation_s.vector + 1) = coordinates[atom_one][1] - coordinates[atom_two][1];
      *(seperation_s.vector + 2) = coordinates[atom_one][2] - coordinates[atom_two][2];
      
      //translate the information in the bond
      
      charge_A = coordinates[atom_one][5];
      charge_B = coordinates[atom_two][5];

      //calculate the interatomic distance

      Ad3DVectorLength(&seperation_s);
      
      //add this length to the linked list wlement the bond belongs to
      
      Interaction->length = seperation_s.length;

/*    if(seperation_s.length > cut)
            return;*/
                  
      //get reciprocal of seperation

      length_rec = 1/seperation_s.length;
      
      //calculate vdw holder
      
      vdw_hold = pow(length_rec, 6);

      //while this is being calculate 
      //calculate some non dependant variables

      //est holder

      est_hold = (EPSILON_RP*charge_A*charge_B*length_rec);

      //est shift factor
      
      shift_fac = (cut - seperation_s.length)*(cut - seperation_s.length)*r_cutoff2;
      
      //vdw holder may be finished now so calculate vdw consts
      
      A = Interaction->params[0]*vdw_hold*vdw_hold;
      B = Interaction->params[1]*vdw_hold;
      
      //modify the potential by the shift function
      
      *est_pot += (est_hold*shift_fac);
      *vdw_pot += A - B;      
}

inline void calc_GRID_EST(ListElement *Interaction, double *est_pot, double EPSILON_RP, double **coordinates,
                        double cut, int grid_point)
{
      int atom_one, atom_two;
      double length_rec, est_hold;
      double charge_A, charge_B;
      Vector3D seperation_s;

      atom_one = Interaction->bond[0];
      atom_two = Interaction->bond[1];
                  
      //calculate seperation vectir (r1 - r2)
      
      *(seperation_s.vector + 0) = coordinates[atom_one][0] - coordinates[atom_two][0];
      *(seperation_s.vector + 1) = coordinates[atom_one][1] - coordinates[atom_two][1];
      *(seperation_s.vector + 2) = coordinates[atom_one][2] - coordinates[atom_two][2];
      
      //translate the information in the bond
      
      charge_A = coordinates[atom_one][5];
      charge_B = coordinates[atom_two][5];

      //calculate the interatomic distance

      Ad3DVectorLength(&seperation_s);
      
      //add this length to the linked list wlement the bond belongs to
      
      Interaction->length = seperation_s.length;

      if(seperation_s.length > cut)
            return;
                  
      //get reciprocal of seperation

      length_rec = 1/seperation_s.length;
      
      //est holder

      est_hold = (EPSILON_RP*charge_A*charge_B*length_rec);

      est_pot[grid_point] += est_hold;
}

inline void grid_EST_switched(ListElement *Interaction, double *est_pot, double EPSILON_RP, double **coordinates,
                        double cutoff_sq, double buffer_sq, int grid_point)
{
      int atom_one, atom_two;
      double length_rec, est_hold, poly;
      double charge_A, charge_B, sep_sq;
      Vector3D seperation_s;

      //could precompute these 
      
      atom_one = Interaction->bond[0];
      atom_two = Interaction->bond[1];
                  
      //calculate seperation vectir (r1 - r2)
      
      *(seperation_s.vector + 0) = coordinates[atom_one][0] - coordinates[atom_two][0];
      *(seperation_s.vector + 1) = coordinates[atom_one][1] - coordinates[atom_two][1];
      *(seperation_s.vector + 2) = coordinates[atom_one][2] - coordinates[atom_two][2];
      
      //translate the information in the bond
      
      charge_A = coordinates[atom_one][5];
      charge_B = coordinates[atom_two][5];

      //calculate the interatomic distance

      Ad3DVectorLength(&seperation_s);
      
      sep_sq = seperation_s.length*seperation_s.length;

      //add this length to the linked list wlement the bond belongs to
      
      Interaction->length = seperation_s.length;

      if(sep_sq > cutoff_sq)
            return;
                  
      //get reciprocal of seperation

      length_rec = 1/seperation_s.length;
      
      //est holder

      est_hold = (EPSILON_RP*charge_A*charge_B*length_rec);

      //if the atom is in the buffer region apply the switch

      if(sep_sq > buffer_sq)
      {
            poly = (sep_sq - buffer_sq)/(cutoff_sq - buffer_sq);
            poly *= poly*(2*poly - 3);
            est_hold *= (1 + poly);
      }
      
      est_pot[grid_point] += est_hold;
}

inline void AdEnzymixGRFNonbondedEnergy(ListElement *Interaction, double *vdw_pot, double *est_pot, double EPSILON_RP, double **coordinates,
                        double cut, double b0)
{
      int atom_one, atom_two;
      double length_rec, vdw_hold, est_hold;
      double A, B;
      double charge_A, charge_B;
      Vector3D seperation_s;
      
      atom_one = Interaction->bond[0];
      atom_two = Interaction->bond[1];
                  
      //calculate seperation vectir (r1 - r2)
      
      *(seperation_s.vector + 0) = coordinates[atom_one][0] - coordinates[atom_two][0];
      *(seperation_s.vector + 1) = coordinates[atom_one][1] - coordinates[atom_two][1];
      *(seperation_s.vector + 2) = coordinates[atom_one][2] - coordinates[atom_two][2];
      
      //translate the information in the bond
      
      charge_A = coordinates[atom_one][5];
      charge_B = coordinates[atom_two][5];

      //calculate the interatomic distance

      Ad3DVectorLength(&seperation_s);
      
      //add this length to the linked list wlement the bond belongs to
      
      Interaction->length = seperation_s.length;

/*    if(seperation_s.length > cut)
            return;*/
                  
      //get reciprocal of seperation

      length_rec = 1/seperation_s.length;
      
      //calculate vdw holder
      
      vdw_hold = pow(length_rec, 6);

      //while this is being calculate 
      //calculate some non dependant variables

      //est holder

      est_hold = EPSILON_RP*charge_A*charge_B*(length_rec + b0);

      //vdw holder may be finished now so calculate vdw consts
      
      A = Interaction->params[0]*vdw_hold*vdw_hold;
      B = Interaction->params[1]*vdw_hold;

/*    if(atom_two==0)
            printf("%-12lf%-12lf%-12lf%-12lf\n", charge_A, charge_B, seperation_s.length, est_hold);*/
      
      *est_pot += est_hold;
      *vdw_pot += A - B;      
}


Generated by  Doxygen 1.6.0   Back to index