]> www.hackdaworld.org Git - physik/posic.git/blob - potentials/harmonic_oscillator.c
added renew feature
[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_set_params(t_moldyn *moldyn,int element) {
23
24         t_ho_params *p;
25
26         // set cutoff before parameters (actually only necessary for some pots)
27         if(moldyn->cutoff==0.0) {
28                 printf("[harmonic oscillator] WARNING: no cutoff!\n");
29                 return -1;
30         }
31
32         /* alloc mem for potential parameters */
33         moldyn->pot_params=malloc(sizeof(t_ho_params));
34         if(moldyn->pot_params==NULL) {
35                 perror("[harmonic oscillator] pot params alloc");
36                 return -1;
37         }
38
39         /* these are now ho parameters */
40         p=moldyn->pot_params;
41
42         switch(element) {
43                 case SI:
44                         /* type: silicon */
45                         p->spring_constant=HO_SC_SI;
46                         p->equilibrium_distance=HO_ED_SI;
47                         break;
48                 case C:
49                         /* type: carbon */
50                         p->spring_constant=HO_SC_C;
51                         p->equilibrium_distance=HO_ED_C;
52                         break;
53                 default:
54                         printf("[harmonic oscillator] WARNING: element\n");
55                         return -1;
56         }
57
58         return 0;
59 }
60
61 int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
62
63         t_ho_params *params;
64         t_3dvec force,distance;
65         double d,f;
66         double sc,equi_dist;
67         double energy;
68
69         params=moldyn->pot_params;
70         sc=params->spring_constant;
71         equi_dist=params->equilibrium_distance;
72
73         if(ai<aj) return 0;
74
75         v3_sub(&distance,&(aj->r),&(ai->r));
76
77         if(bc) check_per_bound(moldyn,&distance);
78         d=v3_norm(&distance);
79         if(d<=moldyn->cutoff) {
80                 energy=(0.5*sc*(d-equi_dist)*(d-equi_dist));
81                 moldyn->energy+=energy;
82                 ai->e+=0.5*energy;
83                 aj->e+=0.5*energy;
84                 /* f = -grad E; grad r_ij = -1 1/r_ij distance */
85                 f=sc*(1.0-equi_dist/d);
86                 v3_scale(&force,&distance,f);
87                 v3_add(&(ai->f),&(ai->f),&force);
88                 virial_calc(ai,&force,&distance);
89                 virial_calc(aj,&force,&distance); /* f and d signe switched */
90                 v3_scale(&force,&distance,-f);
91                 v3_add(&(aj->f),&(aj->f),&force);
92         }
93
94         return 0;
95 }
96
97 int harmonic_oscillator_check_2b_bond(t_moldyn *moldyn,
98                                       t_atom *ai,t_atom *aj,u8 bc) {
99
100         t_3dvec distance;
101         double d;
102
103         v3_sub(&distance,&(aj->r),&(ai->r));
104         if(bc) check_per_bound(moldyn,&distance);
105         d=v3_norm(&distance);
106
107         if(d>moldyn->cutoff)
108                 return FALSE;
109
110         return TRUE;
111 }
112