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(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
  23
  24         t_lj_params *params;
  25         t_3dvec force,distance;
  26         double d,h1,h2;
  27         double eps,sig6,sig12;
  28         double energy;
  29
  30         params=moldyn->pot_params;
  31         eps=params->epsilon4;
  32         sig6=params->sigma6;
  33         sig12=params->sigma12;
  34
  35         if(ai<aj) return 0;
  36
  37         v3_sub(&distance,&(aj->r),&(ai->r));
  38         if(bc) check_per_bound(moldyn,&distance);
  39         d=v3_absolute_square(&distance);        /* r^2 */
  40         if(d<=moldyn->cutoff_square) {
  41                 d=1.0/d;                        /* 1/r^2 */
  42                 h2=d*d;                         /* 1/r^4 */
  43                 h2*=d;                          /* 1/r^6 */
  44                 h1=h2*h2;                       /* 1/r^12 */
  45                 energy=(eps*(sig12*h1-sig6*h2)-params->uc);
  46                 moldyn->energy+=energy;         /* total energy */
  47                 ai->e+=0.5*energy;              /* site energy */
  48                 aj->e+=0.5*energy;
  49                 h2*=d;                          /* 1/r^8 */
  50                 h1*=d;                          /* 1/r^14 */
  51                 h2*=6*sig6;
  52                 h1*=12*sig12;
  53                 d=+h1-h2;
  54                 d*=eps;
  55                 v3_scale(&force,&distance,d);
  56                 v3_add(&(aj->f),&(aj->f),&force);
  57                 v3_scale(&force,&force,-1.0); /* f = - grad E */
  58                 v3_add(&(ai->f),&(ai->f),&force);
  59                 virial_calc(ai,&force,&distance);
  60                 virial_calc(aj,&force,&distance); /* f and d signe switched */
  61         }
  62
  63         return 0;
  64 }
  65