From: hackbard Date: Mon, 3 Apr 2006 15:48:42 +0000 (+0000) Subject: added harmonic potntial + bugfixes + boundings X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=45b27e01673a6cc5bebecb49c51d7f587917483e;p=physik%2Fposic.git added harmonic potntial + bugfixes + boundings --- diff --git a/clean b/clean index 2e1cd0e..34c2b3c 100755 --- a/clean +++ b/clean @@ -1,3 +1 @@ -#!/bin/bash - rm -f saves/* video/* diff --git a/moldyn.c b/moldyn.c index 3d669d2..72a672a 100644 --- a/moldyn.c +++ b/moldyn.c @@ -177,6 +177,19 @@ t_3dvec get_total_p(t_atom *atom, int count) { return p_total; } +double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) { + + double tau; + + tau=0.05*nn_dist/(sqrt(3.0*K_BOLTZMANN*t/moldyn->atom[0].mass)); + tau*=1.0E-9; + if(tautau) + printf("[moldyn] warning: time step (%f > %.15f)\n", + moldyn->tau,tau); + + return tau; +} + /* * @@ -189,6 +202,9 @@ t_3dvec get_total_p(t_atom *atom, int count) { int moldyn_integrate(t_moldyn *moldyn) { int i; + int write; + + write=moldyn->write; /* calculate initial forces */ moldyn->force(moldyn); @@ -198,8 +214,7 @@ int moldyn_integrate(t_moldyn *moldyn) { moldyn->integrate(moldyn); /* check for visualiziation */ - // to be continued ... - if(!(i%1)) { + if(!(i%write)) { visual_atoms(moldyn->visual,i*moldyn->tau, moldyn->atom,moldyn->count); } @@ -233,7 +248,7 @@ int velocity_verlet(t_moldyn *moldyn) { /* velocities */ v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass); - v3_add(&(atom[i].r),&(atom[i].r),&delta); + v3_add(&(atom[i].v),&(atom[i].v),&delta); } /* forces depending on chosen potential */ @@ -255,6 +270,72 @@ int velocity_verlet(t_moldyn *moldyn) { * */ +/* harmonic oscillator potential and force */ + +double potential_harmonic_oscillator(t_moldyn *moldyn) { + + t_ho_params *params; + t_atom *atom; + int i,j; + int count; + t_3dvec distance; + double d,u; + double sc,equi_dist; + + params=moldyn->pot_params; + atom=moldyn->atom; + sc=params->spring_constant; + equi_dist=params->equilibrium_distance; + count=moldyn->count; + + u=0.0; + for(i=0;iatom; + count=moldyn->count; + params=moldyn->pot_params; + sc=params->spring_constant; + equi_dist=params->equilibrium_distance; + + for(i=0;idim)); + d=v3_norm(&distance); + if(d<=moldyn->cutoff) { + v3_scale(&force,&distance, + (-sc*(1.0-(equi_dist/d)))); + v3_add(&(atom[i].f),&(atom[i].f),&force); + v3_sub(&(atom[j].f),&(atom[j].f),&force); + } + } + } + + return 0; +} + + /* lennard jones potential & force for one sort of atoms */ double potential_lennard_jones(t_moldyn *moldyn) { diff --git a/moldyn.h b/moldyn.h index d94e193..a2633ba 100644 --- a/moldyn.h +++ b/moldyn.h @@ -20,6 +20,7 @@ typedef struct s_atom { t_3dvec f; /* forces */ int element; /* number of element in pse */ double mass; /* atom mass */ + //t_list vicinity /* verlet neighbour list */ } t_atom; typedef struct s_moldyn { @@ -28,15 +29,22 @@ typedef struct s_moldyn { double (*potential)(struct s_moldyn *moldyn); void *pot_params; int (*force)(struct s_moldyn *moldyn); + double cutoff; double cutoff_square; t_3dvec dim; int (*integrate)(struct s_moldyn *moldyn); int time_steps; double tau; void *visual; + int write; unsigned char status; } t_moldyn; +typedef struct s_ho_params { + double spring_constant; + double equilibrium_distance; +} t_ho_params; + typedef struct s_lj_params { double sigma6; double sigma12; @@ -54,18 +62,20 @@ typedef struct s_lj_params { /* phsical values */ -//#define K_BOLTZMANN 1.3807E-23 -#define K_BOLTZMANN 1.0 +#define K_BOLTZMANN 1.3807e-27 /* Nm/K */ +#define AMU 1.660540e-27 /* kg */ #define FCC 0x01 #define DIAMOND 0x02 #define C 0x06 -#define M_C 6.0 +#define M_C (12.011*AMU) -#define Si 0x0e -#define M_SI 14.0 -#define LC_SI 5.43105 +#define SI 0x0e +#define LC_SI 0.543105e-9 /* m */ +#define M_SI (28.085*AMU) /* kg */ +#define LJ_SIGMA_SI ((0.25*sqrt(3.0)*LC_SI)/1.122462) /* m */ +#define LJ_EPSILON_SI (2.1678*1.60e-19) /* Nm */ /* function prototypes */ @@ -79,9 +89,13 @@ double get_e_pot(t_moldyn *moldyn); double get_total_energy(t_moldyn *moldyn); t_3dvec get_total_p(t_atom *atom,int count); +double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t); + int moldyn_integrate(t_moldyn *moldyn); int velocity_verlet(t_moldyn *moldyn); +double potential_harmonic_oscillator(t_moldyn *moldyn); +int force_harmonic_oscillator(t_moldyn *moldyn); double potential_lennard_jones(t_moldyn *moldyn); int force_lennard_jones(t_moldyn *moldyn); diff --git a/perms b/perms index d1f5e41..5693501 100755 --- a/perms +++ b/perms @@ -1,3 +1 @@ -#!/bin/bash - chmod 640 saves/* diff --git a/posic.c b/posic.c index 926e20c..5ff3bc1 100644 --- a/posic.c +++ b/posic.c @@ -1,9 +1,11 @@ /* * posic.c - precipitation process of silicon carbide in silicon * - * author: Frank Zirkelbach + * author: Frank Zirkelbach * */ + +#include #include "moldyn.h" #include "math/math.h" @@ -29,8 +31,9 @@ int main(int argc,char **argv) { int count; t_lj_params lj; + t_ho_params ho; - char fb[32]="saves/fcc_test"; + char fb[32]="saves/lj_test"; /* init */ @@ -48,6 +51,10 @@ int main(int argc,char **argv) { b=LEN_Y; c=LEN_Z; + vis.dim.x=a*LC_SI; + vis.dim.y=b*LC_SI; + vis.dim.z=c*LC_SI; + t=TEMPERATURE; printf("placing silicon atoms ... "); @@ -55,16 +62,16 @@ int main(int argc,char **argv) { //printf("(%d) ok!\n",count); count=2; si=malloc(2*sizeof(t_atom)); - si[0].r.x=2.0; + si[0].r.x=0.16e-9; si[0].r.y=0; si[0].r.z=0; - si[0].element=Si; - si[0].mass=14.0; - si[1].r.x=-2.0; + si[0].element=SI; + si[0].mass=M_SI; + si[1].r.x=-0.16e-9; si[1].r.y=0; si[1].r.z=0; - si[1].element=Si; - si[1].mass=14.0; + si[1].element=SI; + si[1].mass=M_SI; printf("setting thermal fluctuations\n"); //thermal_init(si,&random,count,t); @@ -74,56 +81,65 @@ int main(int argc,char **argv) { /* check kinetic energy */ e=get_e_kin(si,count); - printf("kinetic energy: %f\n",e); - printf("3/2 N k T = %f\n",1.5*count*K_BOLTZMANN*t); + printf("kinetic energy: %.40f [J]\n",e); + printf("3/2 N k T = %.40f [J]\n",1.5*count*K_BOLTZMANN*t); /* check total momentum */ p=get_total_p(si,count); - printf("total momentum: %f\n",v3_norm(&p)); + printf("total momentum: %.30f [Ns]\n",v3_norm(&p)); /* check potential energy */ md.count=count; md.atom=si; md.potential=potential_lennard_jones; + //md.potential=potential_harmonic_oscillator; md.force=force_lennard_jones; + //md.force=force_harmonic_oscillator; //md.cutoff_square=((LC_SI/4.0)*(LC_SI/4.0)); - md.cutoff_square=36.0; + md.cutoff=(0.4e-9); + md.cutoff_square=(0.6e-9*0.6e-9); md.pot_params=&lj; + //md.pot_params=&ho; md.integrate=velocity_verlet; md.time_steps=RUNS; md.tau=TAU; md.status=0; md.visual=&vis; + md.write=WRITE_FILE; - lj.sigma6=3.0/16.0*LC_SI*LC_SI; + lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI; help=lj.sigma6*lj.sigma6; lj.sigma6*=help; lj.sigma12=lj.sigma6*lj.sigma6; - lj.epsilon=10000; + lj.epsilon=LJ_EPSILON_SI; + + ho.equilibrium_distance=0.3e-9; + ho.spring_constant=1.0e-9; u=get_e_pot(&md); - printf("potential energy: %f\n",u); - printf("total energy (1): %f\n",e+u); - printf("total energy (2): %f\n",get_total_energy(&md)); + printf("potential energy: %.40f [J]\n",u); + printf("total energy (1): %.40f [J]\n",e+u); + printf("total energy (2): %.40f [J]\n",get_total_energy(&md)); md.dim.x=a*LC_SI; md.dim.y=b*LC_SI; md.dim.z=c*LC_SI; + printf("estimated accurate time step: %.30f [s]\n", + estimate_time_step(&md,3.0,t)); + + /* * let's do the actual md algorithm now * * integration of newtons equations */ - /* visualize */ - //visual_atoms(&vis,0.0,si,count); - - moldyn_integrate(&md); - printf("total energy (after integration): %f\n",get_total_energy(&md)); + printf("total energy (after integration): %.40f [J]\n", + get_total_energy(&md)); /* close */ diff --git a/posic.h b/posic.h index 6a8131d..d31331f 100644 --- a/posic.h +++ b/posic.h @@ -17,15 +17,16 @@ #ifndef POSIC_H #define POSIC_H -#define RUNS 200 -#define TAU 0.001 +#define RUNS 10000000 +#define TAU 0.000000000000001 +#define WRITE_FILE 100000 #define TEMPERATURE 273.0 -#define LEN_X 15 -#define LEN_Y 15 -#define LEN_Z 15 +#define LEN_X 3 +#define LEN_Y 3 +#define LEN_Z 3 -#define R_CUTOFF 20 +#define R_CUTOFF (0.25*LEN_Z) #endif diff --git a/visual/visual.c b/visual/visual.c index 9d0f5a2..38ff0b0 100644 --- a/visual/visual.c +++ b/visual/visual.c @@ -31,6 +31,8 @@ int visual_init(t_visual *v,char *filebase) { return -1; } + memset(&(v->dim),0,sizeof(t_3dvec)); + return 0; } @@ -45,9 +47,14 @@ int visual_atoms(t_visual *v,double time,t_atom *atom,int n) { int i,fd; char file[128+64]; + t_3dvec dim; + + dim.x=10e9*v->dim.x; + dim.y=10e9*v->dim.y; + dim.z=10e9*v->dim.z; char pse[19][4]={ - "foo", + "*", "H", "He", "Li", @@ -77,9 +84,9 @@ int visual_atoms(t_visual *v,double time,t_atom *atom,int n) { /* script file update */ dprintf(v->fd,"load xyz %s\n",file); - dprintf(v->fd,"spacefill 200\n"); - //dprintf(v->fd,"rotate x 100\n"); - //dprintf(v->fd,"rotate y 10\n"); + dprintf(v->fd,"spacefill 300\n"); + dprintf(v->fd,"rotate x 100\n"); + dprintf(v->fd,"rotate y 10\n"); dprintf(v->fd,"set ambient 20\n"); dprintf(v->fd,"set specular on\n"); sprintf(file,"%s-%.15f.ppm",v->fb,time); @@ -87,11 +94,23 @@ int visual_atoms(t_visual *v,double time,t_atom *atom,int n) { dprintf(v->fd,"zap\n"); /* write the actual data file */ - dprintf(fd,"%d\n",n); + dprintf(fd,"%d\n",(dim.x==0)?n:n+8); dprintf(fd,"atoms at time %.15f\n",time); for(i=0;i