2 * lennard_jones.c - lennard jones potential
4 * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
12 #include <sys/types.h>
18 #include "../moldyn.h"
19 #include "../math/math.h"
20 #include "lennard_jones.h"
22 int lennard_jones_set_params(t_moldyn *moldyn,int element) {
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");
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;
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");
44 /* these are now lennard jones parameters */
50 p->sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI*LJ_SIGMA_SI;
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);
61 p->sigma6=LJ_SIGMA_C*LJ_SIGMA_C*LJ_SIGMA_C;
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);
71 printf("[lennard jones] WARNING: element\n");
78 int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
81 t_3dvec force,distance;
83 double eps,sig6,sig12;
86 params=moldyn->pot_params;
89 sig12=params->sigma12;
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) {
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 */
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 */
122 int lennard_jones_check_2b_bond(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
127 v3_sub(&dist,&(aj->r),&(ai->r));
128 if(bc) check_per_bound(moldyn,&dist);
129 d=v3_absolute_square(&dist);
131 if(d>moldyn->cutoff_square)