potentials/lennard_jones.c

   1 /*
   2  * lennard_jones.c - lennard jones potential
   3  *
   4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
   5  *
   6  */
   7
   8 #define _GNU_SOURCE
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <string.h>
  12 #include <sys/types.h>
  13 #include <sys/stat.h>
  14 #include <fcntl.h>
  15 #include <unistd.h>
  16 #include <math.h>
  17
  18 #include "../moldyn.h"
  19 #include "../math/math.h"
  20 #include "lennard_jones.h"
  21
  22 int lennard_jones_set_params(t_moldyn *moldyn,int element) {
  23
  24         t_lj_params *p;
  25         t_atom a,b;
  26
  27         // set cutoff before parameters (actually only necessary for some pots)
  28         if(moldyn->cutoff==0.0) {
  29                 printf("[lennard jones] WARNING: no cutoff!\n");
  30                 return -1;
  31         }
  32
  33         /* atoms for correction term */
  34         a.r.x=0.0; a.r.y=0.0; a.r.z=0.0;
  35         b.r.x=0.0; b.r.y=0.0; b.r.z=moldyn->cutoff;
  36
  37         /* alloc mem for potential parameters */
  38         moldyn->pot_params=malloc(sizeof(t_lj_params));
  39         if(moldyn->pot_params==NULL) {
  40                 perror("[lennard jones] pot params alloc");
  41                 return -1;
  42         }
  43
  44         /* these are now lennard jones parameters */
  45         p=moldyn->pot_params;
  46
  47         switch(element) {
  48                 case SI:
  49                         /* type: silicon */
  50                         p->sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI*LJ_SIGMA_SI;
  51                         p->sigma6+=p->sigma6;
  52                         p->sigma12=p->sigma6*p->sigma6;
  53                         p->epsilon4=4.0*LJ_EPSILON_SI;
  54                         p->uc=0.0; // calc it now!
  55                         lennard_jones(moldyn,&a,&b,0);
  56                         p->uc=moldyn->energy;
  57                         moldyn->energy=0.0;
  58                         break;
  59                 case C:
  60                         /* type: carbon */
  61                         p->sigma6=LJ_SIGMA_C*LJ_SIGMA_C*LJ_SIGMA_C;
  62                         p->sigma6+=p->sigma6;
  63                         p->sigma12=p->sigma6*p->sigma6;
  64                         p->epsilon4=4.0*LJ_EPSILON_C;
  65                         p->uc=0.0; // calc it now!
  66                         lennard_jones(moldyn,&a,&b,0);
  67                         p->uc=moldyn->energy;
  68                         moldyn->energy=0.0;
  69                         break;
  70                 default:
  71                         printf("[lennard jones] WARNING: element\n");
  72                         return -1;
  73         }
  74
  75         return 0;
  76 }
  77
  78 int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
  79
  80         t_lj_params *params;
  81         t_3dvec force,distance;
  82         double d,h1,h2;
  83         double eps,sig6,sig12;
  84         double energy;
  85
  86         params=moldyn->pot_params;
  87         eps=params->epsilon4;
  88         sig6=params->sigma6;
  89         sig12=params->sigma12;
  90
  91         if(ai<aj) return 0;
  92
  93         v3_sub(&distance,&(aj->r),&(ai->r));
  94         if(bc) check_per_bound(moldyn,&distance);
  95         d=v3_absolute_square(&distance);        /* r^2 */
  96         if(d<=moldyn->cutoff_square) {
  97                 d=1.0/d;                        /* 1/r^2 */
  98                 h2=d*d;                         /* 1/r^4 */
  99                 h2*=d;                          /* 1/r^6 */
 100                 h1=h2*h2;                       /* 1/r^12 */
 101                 energy=(eps*(sig12*h1-sig6*h2)-params->uc);
 102                 moldyn->energy+=energy;         /* total energy */
 103                 ai->e+=0.5*energy;              /* site energy */
 104                 aj->e+=0.5*energy;
 105                 h2*=d;                          /* 1/r^8 */
 106                 h1*=d;                          /* 1/r^14 */
 107                 h2*=6*sig6;
 108                 h1*=12*sig12;
 109                 d=+h1-h2;
 110                 d*=eps;
 111                 v3_scale(&force,&distance,d);
 112                 v3_add(&(aj->f),&(aj->f),&force);
 113                 v3_scale(&force,&force,-1.0); /* f = - grad E */
 114                 v3_add(&(ai->f),&(ai->f),&force);
 115                 virial_calc(ai,&force,&distance);
 116                 virial_calc(aj,&force,&distance); /* f and d signe switched */
 117         }
 118
 119         return 0;
 120 }
 121
 122 int lennard_jones_check_2b_bond(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 123
 124          t_3dvec dist;
 125          double d;
 126
 127         v3_sub(&dist,&(aj->r),&(ai->r));
 128         if(bc) check_per_bound(moldyn,&dist);
 129         d=v3_absolute_square(&dist);
 130
 131         if(d>moldyn->cutoff_square)
 132                 return FALSE;
 133
 134         return TRUE;
 135 }