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