]> hackdaworld.org Git - physik/posic.git/blob - potentials/harmonic_oscillator.c
introduce a 2body bond function with callback, modified pair corr calc
[physik/posic.git] / potentials / harmonic_oscillator.c
1 /*
2  * harmonic_oscillator.c - harmonic oscillator 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 "harmonic_oscillator.h"
21
22 int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
23
24         t_ho_params *params;
25         t_3dvec force,distance;
26         double d,f;
27         double sc,equi_dist;
28         double energy;
29
30         params=moldyn->pot_params;
31         sc=params->spring_constant;
32         equi_dist=params->equilibrium_distance;
33
34         if(ai<aj) return 0;
35
36         v3_sub(&distance,&(aj->r),&(ai->r));
37
38         if(bc) check_per_bound(moldyn,&distance);
39         d=v3_norm(&distance);
40         if(d<=moldyn->cutoff) {
41                 energy=(0.5*sc*(d-equi_dist)*(d-equi_dist));
42                 moldyn->energy+=energy;
43                 ai->e+=0.5*energy;
44                 aj->e+=0.5*energy;
45                 /* f = -grad E; grad r_ij = -1 1/r_ij distance */
46                 f=sc*(1.0-equi_dist/d);
47                 v3_scale(&force,&distance,f);
48                 v3_add(&(ai->f),&(ai->f),&force);
49                 virial_calc(ai,&force,&distance);
50                 virial_calc(aj,&force,&distance); /* f and d signe switched */
51                 v3_scale(&force,&distance,-f);
52                 v3_add(&(aj->f),&(aj->f),&force);
53         }
54
55         return 0;
56 }
57