From a114a486ccca019135a8cb229545d0e203b61a5b Mon Sep 17 00:00:00 2001 From: hackbard Date: Tue, 28 Nov 2006 13:17:59 +0000 Subject: [PATCH] still: moldyn.c:1236: error: parse error at end of input - :/ --- moldyn.c | 187 ++++++++++++++++++++++++++++++++----------------------- moldyn.h | 31 +++++---- 2 files changed, 130 insertions(+), 88 deletions(-) diff --git a/moldyn.c b/moldyn.c index c0b99eb..00a9ed9 100644 --- a/moldyn.c +++ b/moldyn.c @@ -51,12 +51,12 @@ int moldyn_shutdown(t_moldyn *moldyn) { int set_int_alg(t_moldyn *moldyn,u8 algo) { - switch(alg) { - case 'MOLDYN_INTEGRATE_VERLET': + switch(algo) { + case MOLDYN_INTEGRATE_VERLET: moldyn->integrate=velocity_verlet; break; default: - printf("unknown integration algorithm: %02x\",alg); + printf("unknown integration algorithm: %02x\n",algo); return -1; } @@ -106,25 +106,26 @@ int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) { return 0; } -int set_potential(t_moldyn *moldyn,u8 type,(int *)(func),void *params) { +int set_potential1b(t_moldyn *moldyn,pf_func1b func,void *params) { - switch(type) { - case MOLDYN_1BP: - moldyn->pf_func1b=func; - moldyn->pot1b_params=params; - break; - case MOLDYN_2BP: - moldyn->pf_func2b=func; - moldyn->pot2b_params=params; - break; - case MOLDYN_3BP: - moldyn->pf_func3b=func; - moldyn->pot3b_params=params; - break; - default: - printf("unknown potential type: %02x\n",type); - return -1; - } + moldyn->func1b=func; + moldyn->pot1b_params=params; + + return 0; +} + +int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params) { + + moldyn->func2b=func; + moldyn->pot2b_params=params; + + return 0; +} + +int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params) { + + moldyn->func3b=func; + moldyn->pot3b_params=params; return 0; } @@ -139,7 +140,7 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,char *fb,int timer) { perror("[moldyn] efd open"); return moldyn->efd; } - dprintf("# moldyn total energy log file\n"); + dprintf(moldyn->efd,"# total energy log file\n"); break; case LOG_TOTAL_MOMENTUM: moldyn->mwrite=timer; @@ -148,7 +149,7 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,char *fb,int timer) { perror("[moldyn] mfd open"); return moldyn->mfd; } - dprintf("# moldyn total momentum log file\n"); + dprintf(moldyn->efd,"# total momentum log file\n"); break; case SAVE_STEP: moldyn->swrite=timer; @@ -353,9 +354,14 @@ double get_e_pot(t_moldyn *moldyn) { return moldyn->energy; } +double update_e_kin(t_moldyn *moldyn) { + + return(get_e_kin(moldyn)); +} + double get_total_energy(t_moldyn *moldyn) { - return(get_e_kin(moldyn)+get_e_pot(moldyn)); + return(moldyn->ekin+moldyn->energy); } t_3dvec get_total_p(t_moldyn *moldyn) { @@ -367,7 +373,7 @@ t_3dvec get_total_p(t_moldyn *moldyn) { atom=moldyn->atom; v3_zero(&p_total); - for(i=0;icount;i++) { v3_scale(&p,&(atom[i].v),atom[i].mass); v3_add(&p_total,&p_total,&p); } @@ -401,9 +407,6 @@ int link_cell_init(t_moldyn *moldyn) { lc=&(moldyn->lc); - /* list log fd */ - lc->listfd=open("/dev/null",O_WRONLY); - /* partitioning the md cell */ lc->nx=moldyn->dim.x/moldyn->cutoff; lc->x=moldyn->dim.x/lc->nx; @@ -418,8 +421,7 @@ int link_cell_init(t_moldyn *moldyn) { printf("initializing linked cells (%d)\n",lc->cells); for(i=0;icells;i++) - //list_init(&(lc->subcell[i]),1); - list_init(&(lc->subcell[i])); + list_init(&(lc->subcell[i]),1); link_cell_update(moldyn); @@ -443,7 +445,7 @@ int link_cell_update(t_moldyn *moldyn) { for(i=0;icells;i++) list_destroy(&(moldyn->lc.subcell[i])); - for(count=0;countcount;count++) { + for(count=0;countcount;count++) { i=(atom[count].r.x+(moldyn->dim.x/2))/lc->x; j=(atom[count].r.y+(moldyn->dim.y/2))/lc->y; k=(atom[count].r.z+(moldyn->dim.z/2))/lc->z; @@ -522,8 +524,6 @@ int link_cell_shutdown(t_moldyn *moldyn) { for(i=0;inx*lc->ny*lc->nz;i++) list_shutdown(&(moldyn->lc.subcell[i])); - if(lc->listfd) close(lc->listfd); - return 0; } @@ -533,7 +533,7 @@ int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) { void *ptr; t_moldyn_schedule *schedule; - schedule=moldyn->schedule; + schedule=&(moldyn->schedule); count=++(schedule->content_count); ptr=realloc(moldyn->schedule.runs,count*sizeof(int)); @@ -572,12 +572,15 @@ int moldyn_set_schedule_hook(t_moldyn *moldyn,void *hook,void *hook_params) { int moldyn_integrate(t_moldyn *moldyn) { int i,sched; - unsigned int e,m,s,d,v; + unsigned int e,m,s,v; t_3dvec p; + t_moldyn_schedule *schedule; int fd; char fb[128]; + schedule=&(moldyn->schedule); + /* initialize linked cell method */ link_cell_init(moldyn); @@ -585,7 +588,6 @@ int moldyn_integrate(t_moldyn *moldyn) { e=moldyn->ewrite; m=moldyn->mwrite; s=moldyn->swrite; - d=moldyn->dwrite; v=moldyn->vwrite; if(!(moldyn->lvstat&MOLDYN_LVSTAT_INITIALIZED)) { @@ -598,14 +600,17 @@ int moldyn_integrate(t_moldyn *moldyn) { moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff; /* calculate initial forces */ - moldyn->potential_force_function(moldyn); + potential_force_calc(moldyn); + + /* zero absolute time */ + moldyn->time=0.0; for(sched=0;schedschedule.content_count;sched++) { - /* setting amont of runs and finite time step size */ + /* setting amount of runs and finite time step size */ moldyn->tau=schedule->tau[sched]; moldyn->tau_square=moldyn->tau*moldyn->tau; - moldyn->timesteps=schedule->runs[sched]; + moldyn->time_steps=schedule->runs[sched]; /* integration according to schedule */ @@ -614,18 +619,23 @@ int moldyn_integrate(t_moldyn *moldyn) { /* integration step */ moldyn->integrate(moldyn); + /* increase absolute time */ + moldyn->time+=moldyn->tau; + /* check for log & visualization */ if(e) { if(!(i%e)) dprintf(moldyn->efd, - "%.15f %.45f\n",i*moldyn->tau, + "%.15f %.45f %.45f %.45f\n", + moldyn->time,update_e_kin(moldyn), + moldyn->energy, get_total_energy(moldyn)); } if(m) { if(!(i%m)) { - p=get_total_p(moldyn->atom,moldyn->count); + p=get_total_p(moldyn); dprintf(moldyn->mfd, - "%.15f %.45f\n",i*moldyn->tau, + "%.15f %.45f\n",moldyn->time, v3_norm(&p)); } } @@ -718,13 +728,12 @@ printf("done\n"); int potential_force_calc(t_moldyn *moldyn) { - int i,count; - t_atom *atom; + int i,j,k,count; + t_atom *atom,*btom,*ktom; t_linkcell *lc; t_list neighbour[27]; - t_list *this; - double u; - u8 bc,bc3; + t_list *this,*thisk,*neighbourk; + u8 bc,bck; int countn,dnlc; count=moldyn->count; @@ -741,7 +750,7 @@ int potential_force_calc(t_moldyn *moldyn) { /* single particle potential/force */ if(atom[i].attr&ATOM_ATTR_1BP) - moldyn->pf_func1b(moldyn,&(atom[i])); + moldyn->func1b(moldyn,&(atom[i])); /* 2 body pair potential/force */ if(atom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) { @@ -773,10 +782,10 @@ int potential_force_calc(t_moldyn *moldyn) { if((btom->attr&ATOM_ATTR_2BP)& (atom[i].attr&ATOM_ATTR_2BP)) - moldyn->pf_func2b(moldyn, - &(atom[i]), - btom, - bc); + moldyn->func2b(moldyn, + &(atom[i]), + btom, + bc); /* 3 body potential/force */ @@ -813,10 +822,12 @@ int potential_force_calc(t_moldyn *moldyn) { if(ktom==&(atom[i])) continue; - moldyn->pf_func3b(moldyn,&(atom[i]),btom,ktom,bck); + moldyn->func3b(moldyn,&(atom[i]),btom,ktom,bck); } while(list_next(thisk)!=\ L_NO_NEXT_ELEMENT); + + } } while(list_next(this)!=L_NO_NEXT_ELEMENT); } @@ -833,20 +844,26 @@ int potential_force_calc(t_moldyn *moldyn) { int check_per_bound(t_moldyn *moldyn,t_3dvec *a) { double x,y,z; + t_3dvec *dim; + + dim=&(moldyn->dim); x=0.5*dim->x; y=0.5*dim->y; z=0.5*dim->z; - if(moldyn->MOLDYN_ATTR_PBX) + if(moldyn->status&MOLDYN_STAT_PBX) { if(a->x>=x) a->x-=dim->x; else if(-a->x>x) a->x+=dim->x; - if(moldyn->MOLDYN_ATTR_PBY) + } + if(moldyn->status&MOLDYN_STAT_PBY) { if(a->y>=y) a->y-=dim->y; else if(-a->y>y) a->y+=dim->y; - if(moldyn->MOLDYN_ATTR_PBZ) + } + if(moldyn->status&MOLDYN_STAT_PBZ) { if(a->z>=z) a->z-=dim->z; else if(-a->z>z) a->z+=dim->z; + } return 0; } @@ -858,7 +875,7 @@ int check_per_bound(t_moldyn *moldyn,t_3dvec *a) { /* harmonic oscillator potential and force */ -int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc)) { +int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_ho_params *params; t_3dvec force,distance; @@ -869,7 +886,7 @@ int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc)) { sc=params->spring_constant; equi_dist=params->equilibrium_distance; - v3_sub(&distance,&(ai->r),&(aj->r); + v3_sub(&distance,&(ai->r),&(aj->r)); v3_per_bound(&distance,&(moldyn->dim)); if(bc) check_per_bound(moldyn,&distance); @@ -890,10 +907,10 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_lj_params *params; t_3dvec force,distance; - double d,h1,h2,u; + double d,h1,h2; double eps,sig6,sig12; - params=moldyn->pot_params; + params=moldyn->pot2b_params; eps=params->epsilon4; sig6=params->sigma6; sig12=params->sigma12; @@ -960,12 +977,17 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_tersoff_mult_params *params; t_tersoff_exchange *exchange; - t_3dvec dist_ij; + t_3dvec dist_ij,force; double d_ij; - double A,B,R,S,lambda; + double A,B,R,S,lambda,mu; + double f_r,df_r; + double f_c,df_c; int num; + double s_r; + double arg; + double scale; - params=moldyn->pot_params; + params=moldyn->pot2b_params; num=ai->bnum; exchange=&(params->exchange); @@ -992,7 +1014,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { lambda=params->lambda[num]; /* more constants depending of atoms i and j, needed in 3bp */ params->exchange.B=&(params->B[num]); - params->exchange.mu=params->mu[num]; + params->exchange.mu=&(params->mu[num]); + mu=params->mu[num]; params->exchange.chi=1.0; } else { @@ -1003,6 +1026,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* more constants depending of atoms i and j, needed in 3bp */ params->exchange.B=&(params->Bmixed); params->exchange.mu=&(params->mu_m); + mu=params->mu_m; params->exchange.chi=params->chi; } @@ -1014,7 +1038,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { if(d_ij>S) return 0; - f_r=A*exp(-lamda*d_ij); + f_r=A*exp(-lambda*d_ij); df_r=-lambda*f_r/d_ij; /* f_a, df_a calc + save for 3bp use */ @@ -1029,9 +1053,9 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { } else { s_r=S-R; - arg=PI*(d_ij-R)/s_r; + arg=M_PI*(d_ij-R)/s_r; f_c=0.5+0.5*cos(arg); - df_c=-0.5*sin(arg)*(PI/(s_r*d_ij)); + df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij)); scale=df_c*f_r+df_r*f_c; v3_scale(&force,&dist_ij,scale); } @@ -1060,15 +1084,24 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_3dvec dist_ij,dist_ik,dist_jk; t_3dvec temp,force; double R,S,s_r; - double d_ij,d_ik,d_jk; + double d_ij,d_ij2,d_ik,d_jk; double f_c,df_c,b_ij,f_a,df_a; - double n,c,d,h,neta,betan,betan_1; + double f_c_ik,df_c_ik,arg; + double scale; + double chi; + double n,c,d,h,beta,betan; + double c2,d2,c2d2; + double numer,denom; double theta,cos_theta,sin_theta; + double d_theta,d_theta1,d_theta2; + double h_cos,h_cos2,d2_h_cos2; + double frac1,bracket1,bracket2,bracket2_n_1,bracket2_n; + double bracket3,bracket3_pow_1,bracket3_pow; int num; - params=moldyn->pot_params; + params=moldyn->pot3b_params; num=ai->bnum; - exchange=params->exchange; + exchange=&(params->exchange); if(!(exchange->run3bp)) return 0; @@ -1098,7 +1131,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { */ - v3_sub(&dist_ik,&(aj->i),&(ak->r)); + v3_sub(&dist_ik,&(ai->r),&(ak->r)); if(bc) check_per_bound(moldyn,&dist_ik); d_ik=v3_norm(&dist_ik); @@ -1123,9 +1156,9 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { } else { s_r=S-R; - arg=PI*(d_ik-R)/s_r; + arg=M_PI*(d_ik-R)/s_r; f_c_ik=0.5+0.5*cos(arg); - df_c_ik=-0.5*sin(arg)*(PI/(s_r*d_ik)); + df_c_ik=-0.5*sin(arg)*(M_PI/(s_r*d_ik)); } v3_sub(&dist_jk,&(aj->r),&(ak->r)); @@ -1146,7 +1179,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { denom=2*d_ij*d_ik; cos_theta=numer/denom; sin_theta=sqrt(1.0-(cos_theta*cos_theta)); - theta=arccos(cos_theta); + theta=acos(cos_theta); d_theta=(-1.0/sqrt(1.0-cos_theta*cos_theta))/(denom*denom); d_theta1=2*denom-numer*2*d_ik/d_ij; d_theta2=2*denom-numer*2*d_ij/d_ik; @@ -1171,12 +1204,12 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { b_ij=chi*bracket3_pow; /* derivation of theta */ - v3_scale(&force,&dist_ij,d1_theta); + v3_scale(&force,&dist_ij,d_theta1); v3_scale(&temp,&dist_ik,d_theta2); v3_add(&force,&force,&temp); /* part 1 of derivation of b_ij */ - v3_scale(&force,sin_theta*2*h_cos*f_c_ik*frac1); + v3_scale(&force,&force,sin_theta*2*h_cos*f_c_ik*frac1); /* part 2 of derivation of b_ij */ v3_scale(&temp,&dist_ik,df_c_ik*bracket1); diff --git a/moldyn.h b/moldyn.h index dafd86c..ae133f7 100644 --- a/moldyn.h +++ b/moldyn.h @@ -70,11 +70,12 @@ typedef struct s_moldyn { t_3dvec dim; /* dimensions of the simulation volume */ /* potential force function and parameter pointers */ - int (*pf_func1b)(struct s_moldyn *,t_atom *); + int (*func1b)(struct s_moldyn *moldyn,t_atom *ai); void *pot1b_params; - int (*pf_func2b)(struct s_moldyn *,t_atom *,t_atom *); + int (*func2b)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); void *pot2b_params; - int (*pf_func3b)(struct s_moldyn *,t_atom *,t_atom *,t_atom *); + int (*func3b)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak, + u8 bck); void *pot3b_params; //int (*potential_force_function)(struct s_moldyn *moldyn); @@ -93,6 +94,7 @@ typedef struct s_moldyn { int (*integrate)(struct s_moldyn *moldyn); int time_steps; /* amount of iterations */ double tau; /* delta t */ + double time; /* absolute time */ double tau_square; /* delta t squared */ double elapsed; /* total elapsed time */ @@ -163,6 +165,7 @@ typedef struct s_lj_params { /* tersoff exchange structure to exchange 2bp and 3bp calculated values */ typedef struct s_tersoff_exchange { double f_c,df_c; + double f_a,df_a; t_3dvec dist_ij; double d_ij; @@ -211,7 +214,7 @@ typedef struct s_tersoff_mult_params { double h[2]; t_tersoff_exchange exchange; /* exchange between 2bp and 3bp calc */ -} t_tersoff_params; +} t_tersoff_mult_params; @@ -235,10 +238,10 @@ typedef struct s_tersoff_mult_params { #define MOLDYN_POTENTIAL_LJ 0x01 #define MOLDYN_POTENTIAL_TM 0x02 -#define MOLDYN_SET_POTENTIAL 0x00 -#define MOLDYN_SET_TEMPERATURE 0x01 -#define MOLDYN_SET_ -#define MOLDYN_SET_ +#define LOG_TOTAL_ENERGY 0x01 +#define LOG_TOTAL_MOMENTUM 0x02 +#define SAVE_STEP 0x04 +#define VISUAL_STEP 0x08 #define TRUE 1 #define FALSE 0 @@ -271,6 +274,10 @@ typedef struct s_tersoff_mult_params { * */ +typedef int (*pf_func1b)(t_moldyn *,t_atom *ai); +typedef int (*pf_func2b)(t_moldyn *,t_atom *,t_atom *,u8 bc); +typedef int (*pf_func3b)(t_moldyn *,t_atom *,t_atom *,t_atom *,u8 bc); + int moldyn_init(t_moldyn *moldyn,int argc,char **argv); int moldyn_shutdown(t_moldyn *moldyn); @@ -279,10 +286,12 @@ int set_cutoff(t_moldyn *moldyn,double cutoff); int set_temperature(t_moldyn *moldyn,double t); int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize); int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z); -int set_potential(t_moldyn *moldyn,u8 type,(int *)(func),void *params); +int set_potential1b(t_moldyn *moldyn,pf_func1b func,void *params); +int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params); +int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params); int moldyn_set_log(t_moldyn *moldyn,u8 type,char *fb,int timer); -nt moldyn_log_shutdown(t_moldyn *moldyn); +int moldyn_log_shutdown(t_moldyn *moldyn); int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, u8 attr,u8 bnum,int a,int b,int c); @@ -313,7 +322,7 @@ int velocity_verlet(t_moldyn *moldyn); int potential_force_calc(t_moldyn *moldyn); int check_per_bound(t_moldyn *moldyn,t_3dvec *a); -int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_Atom *aj,u8 bc); +int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai); int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); -- 2.20.1