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