2 * harmonic_oscillator.c - harmonic oscillator 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 "harmonic_oscillator.h"
22 int harmonic_oscillator_set_params(t_moldyn *moldyn,int element) {
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");
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");
39 /* these are now ho parameters */
45 p->spring_constant=HO_SC_SI;
46 p->equilibrium_distance=HO_ED_SI;
50 p->spring_constant=HO_SC_C;
51 p->equilibrium_distance=HO_ED_C;
54 printf("[harmonic oscillator] WARNING: element\n");
62 int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
65 t_3dvec force,distance;
70 params=moldyn->pot_params;
71 sc=params->spring_constant;
72 equi_dist=params->equilibrium_distance;
76 v3_sub(&distance,&(aj->r),&(ai->r));
78 if(bc) check_per_bound(moldyn,&distance);
80 if(d<=moldyn->cutoff) {
81 energy=(0.5*sc*(d-equi_dist)*(d-equi_dist));
82 moldyn->energy+=energy;
85 /* f = -grad E; grad r_ij = -1 1/r_ij distance */
86 f=sc*(1.0-equi_dist/d);
87 v3_scale(&force,&distance,f);
88 v3_add(&(ai->f),&(ai->f),&force);
89 virial_calc(ai,&force,&distance);
90 virial_calc(aj,&force,&distance); /* f and d signe switched */
91 v3_scale(&force,&distance,-f);
92 v3_add(&(aj->f),&(aj->f),&force);
98 int harmonic_oscillator_check_2b_bond(t_moldyn *moldyn,
99 t_atom *ai,t_atom *aj,u8 bc) {
104 v3_sub(&distance,&(aj->r),&(ai->r));
105 if(bc) check_per_bound(moldyn,&distance);
106 d=v3_norm(&distance);