]> hackdaworld.org Git - physik/posic.git/blob - potentials/lennard_jones.c
some fixes
[physik/posic.git] / 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);        /* r^2 */
39         if(d<=moldyn->cutoff_square) {
40                 d=1.0/d;                        /* 1/r^2 */
41                 h2=d*d;                         /* 1/r^4 */
42                 h2*=d;                          /* 1/r^6 */
43                 h1=h2*h2;                       /* 1/r^12 */
44                 moldyn->energy+=(eps*(sig12*h1-sig6*h2)-params->uc);
45                 h2*=d;                          /* 1/r^8 */
46                 h1*=d;                          /* 1/r^14 */
47                 h2*=6*sig6;
48                 h1*=12*sig12;
49                 d=+h1-h2;
50                 d*=eps;
51                 v3_scale(&force,&distance,d);
52                 v3_add(&(aj->f),&(aj->f),&force);
53                 v3_scale(&force,&force,-1.0); /* f = - grad E */
54                 v3_add(&(ai->f),&(ai->f),&force);
55                 virial_calc(ai,&force,&distance);
56                 virial_calc(aj,&force,&distance); /* f and d signe switched */
57         }
58
59         return 0;
60 }
61