From 792f14f95b47989f7f12df0ea70b54619be016ee Mon Sep 17 00:00:00 2001 From: hackbard Date: Tue, 4 Apr 2006 10:09:32 +0000 Subject: [PATCH] checkin before cafeteria visit ;) --- moldyn.c | 200 +++++++++++++++++++++++++++++++++++++++++++++++++++---- moldyn.h | 45 +++++++++++-- posic.c | 159 ++++--------------------------------------- 3 files changed, 240 insertions(+), 164 deletions(-) diff --git a/moldyn.c b/moldyn.c index aec87f6..f2ac637 100644 --- a/moldyn.c +++ b/moldyn.c @@ -5,17 +5,136 @@ * */ -#include "moldyn.h" - +#define _GNU_SOURCE #include #include +#include +#include +#include +#include +#include #include +#include "moldyn.h" + #include "math/math.h" #include "init/init.h" #include "random/random.h" #include "visual/visual.h" +int moldyn_usage(char **argv) { + + printf("\n%s usage:\n\n",argv[0]); + printf("--- general options ---\n"); + printf("-E (log total energy)\n"); + printf("-M (log total momentum)\n"); + printf("-D (dump total information)\n"); + printf("-S (single save file)\n"); + printf("--- physics options ---\n"); + printf("-T [K] (%f)\n",MOLDYN_TEMP); + printf("-t [s] (%f)\n",MOLDYN_TAU); + printf("-R (%d)\n",MOLDYN_RUNS); + printf("\n"); + + return 0; +} + +int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) { + + int i; + + memset(moldyn,0,sizeof(t_moldyn)); + + /* default values */ + moldyn->t=MOLDYN_TEMP; + moldyn->tau=MOLDYN_TAU; + moldyn->time_steps=MOLDYN_RUNS; + + /* parse argv */ + for(i=1;iewrite=atoi(argv[++i]); + strncpy(moldyn->efb,argv[++i],64); + break; + case 'M': + moldyn->mwrite=atoi(argv[++i]); + strncpy(moldyn->mfb,argv[++i],64); + break; + case 'D': + moldyn->dwrite=atoi(argv[++i]); + strncpy(moldyn->dfb,argv[++i],64); + break; + case 'S': + moldyn->swrite=atoi(argv[++i]); + strncpy(moldyn->sfb,argv[++i],64); + break; + case 'T': + break; + moldyn->t=atof(argv[++i]); + case 't': + moldyn->tau=atof(argv[++i]); + break; + case 'R': + moldyn->time_steps=atoi(argv[++i]); + break; + default: + printf("unknown option %s\n",argv[i]); + moldyn_usage(argv); + return -1; + } + } else { + moldyn_usage(argv); + return -1; + } + } + + return 0; +} + +int moldyn_log_init(t_moldyn *moldyn) { + + moldyn->lvstat=0; + + if(moldyn->ewrite) { + moldyn->efd=open(moldyn->efb,O_WRONLY|O_CREAT|O_TRUNC); + if(moldyn->efd<0) { + perror("[moldyn] efd open"); + return moldyn->efd; + } + dprintf(moldyn->efd,"# moldyn total energy logfile\n"); + moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_E; + } + + if(moldyn->mwrite) { + moldyn->mfd=open(moldyn->mfb,O_WRONLY|O_CREAT|O_TRUNC); + if(moldyn->mfd<0) { + perror("[moldyn] mfd open"); + return moldyn->mfd; + } + dprintf(moldyn->mfd,"# moldyn total momentum logfile\n"); + moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_M; + } + + if(moldyn->swrite) + moldyn->lvstat|=MOLDYN_LVSTAT_SAVE; + + if(moldyn->dwrite) { + moldyn->dfd=open(moldyn->dfb,O_WRONLY|O_CREAT|O_TRUNC); + if(moldyn->dfd<0) { + perror("[moldyn] dfd open"); + return moldyn->dfd; + } + write(moldyn->dfd,moldyn,sizeof(t_moldyn)); + moldyn->lvstat|=MOLDYN_LVSTAT_DUMP; + } + + if(moldyn->dwrite) + moldyn->lvstat|=MOLDYN_LVSTAT_VISUAL; + + return 0; +} int create_lattice(unsigned char type,int element,double mass,double lc, int a,int b,int c,t_atom **atom) { @@ -73,7 +192,7 @@ int destroy_lattice(t_atom *atom) { return 0; } -int thermal_init(t_atom *atom,t_random *random,int count,double t) { +int thermal_init(t_moldyn *moldyn,t_random *random,int count) { /* * - gaussian distribution of velocities @@ -84,11 +203,14 @@ int thermal_init(t_atom *atom,t_random *random,int count,double t) { int i; double v,sigma; t_3dvec p_total,delta; + t_atom *atom; + + atom=moldyn->atom; /* gaussian distribution of velocities */ v3_zero(&p_total); for(i=0;it/atom[i].mass); /* x direction */ v=sigma*rand_get_gauss(random); atom[i].v.x=v; @@ -111,15 +233,18 @@ int thermal_init(t_atom *atom,t_random *random,int count,double t) { } /* velocity scaling */ - scale_velocity(atom,count,t); + scale_velocity(moldyn,count); return 0; } -int scale_velocity(t_atom *atom,int count,double t) { +int scale_velocity(t_moldyn *moldyn,int count) { int i; double e,c; + t_atom *atom; + + atom=moldyn->atom; /* * - velocity scaling (E = 3/2 N k T), E: kinetic energy @@ -127,7 +252,7 @@ int scale_velocity(t_atom *atom,int count,double t) { e=0.0; for(i=0;it)); for(i=0;iwrite; + e=moldyn->ewrite; + m=moldyn->mwrite; + s=moldyn->swrite; + d=moldyn->dwrite; + v=moldyn->vwrite; + + if(!(lvstat&MOLDYN_LVSTAT_INITIALIZED)) { + printf("[moldyn] warning, lv system not initialized\n"); + return -1; + } /* calculate initial forces */ moldyn->force(moldyn); @@ -213,11 +352,44 @@ int moldyn_integrate(t_moldyn *moldyn) { /* integration step */ moldyn->integrate(moldyn); - /* check for visualiziation */ - if(!(i%write)) { - visual_atoms(moldyn->visual,i*moldyn->tau, - moldyn->atom,moldyn->count); - printf("finished %d / %d\n",i,moldyn->time_steps); + /* check for log & visualiziation */ + if(e) { + if(!(i%e)) + dprintf(moldyn->efd, + "%.15f %.45f\n",i*moldyn->tau, + get_total_energy(moldyn)); + } + if(m) { + if(!(i%m)) { + p=get_total_p(moldyn->atom,moldyn->count); + dprintf(moldyn->mfd, + "%.15f %.45f\n",i*moldyn->tau, + v3_norm(&p)); + } + } + if(s) { + if(!(i%s)) { + snprintf(fb,128,"%s-%f-%.15f.save",moldyn->sfb, + moldyn->t,i*moldyn->tau); + fd=open(fb,O_WRONLY|O_TRUNC|O_CREAT); + if(fd<0) perror("[moldyn] save fd open"); + else { + write(fd,moldyn,sizeof(t_moldyn)); + write(fd,moldyn->atom, + moldyn->count*sizeof(t_atom)); + } + } + } + if(d) { + if(!(i%d)) + write(moldyn->dfd,moldyn->atom, + moldyn->count*sizeof(t_atom)); + + } + if(v) { + if(!(i%v)) + visual_atoms(moldyn->visual,i*moldyn->tau, + moldyn->atom,moldyn->count); } } diff --git a/moldyn.h b/moldyn.h index a2633ba..cbfc173 100644 --- a/moldyn.h +++ b/moldyn.h @@ -24,19 +24,39 @@ typedef struct s_atom { } t_atom; typedef struct s_moldyn { + /* atoms, amount, dimensions */ int count; t_atom *atom; + t_3dvec dim; + /* potential, force & parameters */ double (*potential)(struct s_moldyn *moldyn); - void *pot_params; int (*force)(struct s_moldyn *moldyn); + void *pot_params; double cutoff; double cutoff_square; - t_3dvec dim; + /* temperature */ + double t; + /* integration of newtons equations */ int (*integrate)(struct s_moldyn *moldyn); int time_steps; double tau; + /* logging & visualization */ + unsigned char lvstat; + unsigned int ewrite; + char efb[64]; + int efd; + unsigned int mwrite; + char mfb[64]; + int mfd; + unsigned int swrite; + char sfb[64]; + int sfd; + unsigned int dwrite; + char dfb[64]; + int dfd; + unsigned int vwrite; void *visual; - int write; + /* moldyn general status */ unsigned char status; } t_moldyn; @@ -57,9 +77,20 @@ typedef struct s_lj_params { /* general defines */ +#define MOLDYN_LVSTAT_TOTAL_E 0x01 +#define MOLDYN_LVSTAT_TOTAL_M 0x02 +#define MOLDYN_LVSTAT_SAVE 0x04 +#define MOLDYN_LVSTAT_DUMP 0x08 +#define MOLDYN_LVSTAT_VISUAL 0x10 +#define MOLDYN_LVSTAT_INITIALIZED 0x80 + #define MOLDYN_STAT_POTENTIAL 0x01 #define MOLDYN_STAT_FORCE 0x02 +#define MOLDYN_TEMP 273.0 +#define MOLDYN_TAU 1.0e-15 +#define MOLDYN_RUNS 1000000 + /* phsical values */ #define K_BOLTZMANN 1.3807e-27 /* Nm/K */ @@ -79,11 +110,15 @@ typedef struct s_lj_params { /* function prototypes */ +int moldyn_usage(char **argv); +int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv); +int moldyn_log_init(t_moldyn *moldyn); + int create_lattice(unsigned char type,int element,double mass,double lc, int a,int b,int c,t_atom **atom); int destroy_lattice(t_atom *atom); -int thermal_init(t_atom *atom,t_random *random,int count,double t); -int scale_velocity(t_atom *atom,int count,double t); +int thermal_init(t_moldyn *moldyn,t_random *random,int count); +int scale_velocity(t_moldyn *moldyn,int count); double get_e_kin(t_atom *atom,int count); double get_e_pot(t_moldyn *moldyn); double get_total_energy(t_moldyn *moldyn); diff --git a/posic.c b/posic.c index 5ff3bc1..7c1c418 100644 --- a/posic.c +++ b/posic.c @@ -31,7 +31,6 @@ int main(int argc,char **argv) { int count; t_lj_params lj; - t_ho_params ho; char fb[32]="saves/lj_test"; @@ -51,15 +50,18 @@ 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; + /* set for 'bounding atoms' */ + //vis.dim.x=a*LC_SI; + //vis.dim.y=b*LC_SI; + //vis.dim.z=c*LC_SI; t=TEMPERATURE; printf("placing silicon atoms ... "); - //count=create_lattice(DIAMOND,Si,M_SI,LC_SI,a,b,c,&si); - //printf("(%d) ok!\n",count); + count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&si); + printf("(%d) ok!\n",count); + + /* testing purpose count=2; si=malloc(2*sizeof(t_atom)); si[0].r.x=0.16e-9; @@ -72,11 +74,12 @@ int main(int argc,char **argv) { si[1].r.z=0; si[1].element=SI; si[1].mass=M_SI; + */ printf("setting thermal fluctuations\n"); - //thermal_init(si,&random,count,t); - v3_zero(&(si[0].v)); - v3_zero(&(si[1].v)); + thermal_init(si,&random,count,t); + //v3_zero(&(si[0].v)); + //v3_zero(&(si[1].v)); /* check kinetic energy */ @@ -92,14 +95,10 @@ int main(int argc,char **argv) { 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=(0.4e-9); - md.cutoff_square=(0.6e-9*0.6e-9); + md.cutoff=R_CUTOFF; + md.cutoff_square=(R_CUTOFF*R_CUTOFF); md.pot_params=&lj; - //md.pot_params=&ho; md.integrate=velocity_verlet; md.time_steps=RUNS; md.tau=TAU; @@ -113,9 +112,6 @@ int main(int argc,char **argv) { lj.sigma12=lj.sigma6*lj.sigma6; 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: %.40f [J]\n",u); @@ -147,133 +143,6 @@ int main(int argc,char **argv) { rand_close(&random); - - //printf("starting velocity verlet: "); - //fflush(stdout); - - //for(runs=0;runsLX) si[i].x-=LEN_X; -// else if(si[i].x<-LX) si[i].x+=LEN_X; -// si[i].y+=(tau2*si[i].fy/m2); -// if(si[i].y>LY) si[i].y-=LEN_Y; -// else if(si[i].y<-LY) si[i].y+=LEN_Y; -// si[i].z+=(tau2*si[i].fz/m2); -// if(si[i].z>LZ) si[i].z-=LEN_Z; -// else if(si[i].z<-LZ) si[i].z+=LEN_Z; -// /* calculation of velocities v(t+h/2) */ -// si[i].vx+=(tau*si[i].fx/m2); -// si[i].vy+=(tau*si[i].fy/m2); -// si[i].vz+=(tau*si[i].fz/m2); -// /* reset of forces */ -// si[i].fx=.0; -// si[i].fy=.0; -// si[i].fz=.0; -// } -// for(i=0;iLX) deltax-=LEN_X; -// else if(-deltax>LX) deltax+=LEN_X; -// deltax2=deltax*deltax; -// deltay=si[i].y-si[j].y; -// if(deltay>LY) deltay-=LEN_Y; -// else if(-deltay>LY) deltay+=LEN_Y; -// deltay2=deltay*deltay; -// deltaz=si[i].z-si[j].z; -// if(deltaz>LZ) deltaz-=LEN_Z; -// else if(-deltaz>LZ) deltaz+=LEN_Z; -// deltaz2=deltaz*deltaz; -// distance=deltax2+deltay2+deltaz2; -// if(distance<=R2_CUTOFF) { -// tmp=1.0/distance; // 1/r^2 -// lj1=tmp; // 1/r^2 -// tmp*=tmp; // 1/r^4 -// lj1*=tmp; // 1/r^6 -// tmp*=tmp; // 1/r^8 -// lj2=tmp; // 1/r^8 -// lj1*=tmp; // 1/r^14 -// lj1*=LJ_SIGMA_12; -// lj2*=LJ_SIGMA_06; -// lj=-2*lj1+lj2; -// si[i].fx-=lj*deltax; -// si[i].fy-=lj*deltay; -// si[i].fz-=lj*deltaz; -// si[j].fx+=lj*deltax; -// si[j].fy+=lj*deltay; -// si[j].fz+=lj*deltaz; -// } -// } -// } -// for(i=0;i