From c0a8b254109929fba10795e644187c51742108a8 Mon Sep 17 00:00:00 2001 From: hackbard Date: Fri, 15 Dec 2006 15:49:06 +0000 Subject: [PATCH] pot hopefully ok, doing timestep tests for 450C now ... --- moldyn.c | 65 +++++++++++++++++++++++++++++++++++++++++++------------- moldyn.h | 3 ++- run | 2 +- sic.c | 30 ++++++++++++++------------ 4 files changed, 70 insertions(+), 30 deletions(-) diff --git a/moldyn.c b/moldyn.c index f8bfd3f..5aa6d87 100644 --- a/moldyn.c +++ b/moldyn.c @@ -767,23 +767,23 @@ int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) { t_moldyn_schedule *schedule; schedule=&(moldyn->schedule); - count=++(schedule->content_count); + count=++(schedule->total_sched); - ptr=realloc(moldyn->schedule.runs,count*sizeof(int)); + ptr=realloc(schedule->runs,count*sizeof(int)); if(!ptr) { perror("[moldyn] realloc (runs)"); return -1; } - moldyn->schedule.runs=ptr; - moldyn->schedule.runs[count-1]=runs; + schedule->runs=ptr; + schedule->runs[count-1]=runs; ptr=realloc(schedule->tau,count*sizeof(double)); if(!ptr) { perror("[moldyn] realloc (tau)"); return -1; } - moldyn->schedule.tau=ptr; - moldyn->schedule.tau[count-1]=tau; + schedule->tau=ptr; + schedule->tau[count-1]=tau; return 0; } @@ -806,16 +806,16 @@ int moldyn_set_schedule_hook(t_moldyn *moldyn,void *hook,void *hook_params) { int moldyn_integrate(t_moldyn *moldyn) { - int i,sched; + int i; unsigned int e,m,s,v; t_3dvec p; - t_moldyn_schedule *schedule; + t_moldyn_schedule *sched; t_atom *atom; int fd; char dir[128]; double ds; - schedule=&(moldyn->schedule); + sched=&(moldyn->schedule); atom=moldyn->atom; /* initialize linked cell method */ @@ -852,12 +852,12 @@ int moldyn_integrate(t_moldyn *moldyn) { moldyn->debug=0; /* executing the schedule */ - for(sched=0;schedschedule.content_count;sched++) { + for(sched->count=0;sched->counttotal_sched;sched->count++) { /* setting amount of runs and finite time step size */ - moldyn->tau=schedule->tau[sched]; + moldyn->tau=sched->tau[sched->count]; moldyn->tau_square=moldyn->tau*moldyn->tau; - moldyn->time_steps=schedule->runs[sched]; + moldyn->time_steps=sched->runs[sched->count]; /* integration according to schedule */ @@ -907,7 +907,7 @@ int moldyn_integrate(t_moldyn *moldyn) { visual_atoms(&(moldyn->vis),moldyn->time, moldyn->atom,moldyn->count); printf("\rsched: %d, steps: %d, debug: %d", - sched,i,moldyn->debug); + sched->count,i,moldyn->debug); fflush(stdout); } } @@ -918,8 +918,8 @@ int moldyn_integrate(t_moldyn *moldyn) { } /* check for hooks */ - if(schedule->hook) - schedule->hook(moldyn,schedule->hook_params); + if(sched->hook) + sched->hook(moldyn,sched->hook_params); /* get a new info line */ printf("\n"); @@ -1108,6 +1108,9 @@ int potential_force_calc(t_moldyn *moldyn) { } } +#ifdef DEBUG +printf("\n\n"); +#endif return 0; } @@ -1397,6 +1400,14 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* add forces of 2bp (ij, ji) contribution * dVij = dVji and we sum up both: no 1/2) */ v3_add(&(ai->f),&(ai->f),&force); +#ifdef DEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVij, dVji (2bp) contrib:\n"); + printf("%f | %f\n",force.x,ai->f.x); + printf("%f | %f\n",force.y,ai->f.y); + printf("%f | %f\n",force.z,ai->f.z); +} +#endif /* energy 2bp contribution (ij, ji) is 0.5 f_r f_c ... */ moldyn->energy+=(0.5*f_r*f_c); @@ -1488,6 +1499,14 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* add force */ v3_add(&(ai->f),&(ai->f),&force); +#ifdef DEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVij (3bp) contrib:\n"); + printf("%f | %f\n",force.x,ai->f.x); + printf("%f | %f\n",force.y,ai->f.y); + printf("%f | %f\n",force.z,ai->f.z); +} +#endif /* add energy of 3bp sum */ moldyn->energy+=(0.5*f_c*b*f_a); @@ -1516,6 +1535,14 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* add force */ v3_add(&(ai->f),&(ai->f),&force); +#ifdef DEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVji (3bp) contrib:\n"); + printf("%f | %f\n",force.x,ai->f.x); + printf("%f | %f\n",force.y,ai->f.y); + printf("%f | %f\n",force.z,ai->f.z); +} +#endif return 0; } @@ -1777,6 +1804,14 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5); v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */ /* scaled with 0.5 ^ */ +#ifdef DEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVjk (3bp) contrib:\n"); + printf("%f | %f\n",temp2.x,ai->f.x); + printf("%f | %f\n",temp2.y,ai->f.y); + printf("%f | %f\n",temp2.z,ai->f.z); +} +#endif } diff --git a/moldyn.h b/moldyn.h index 9493f5f..38e4999 100644 --- a/moldyn.h +++ b/moldyn.h @@ -57,7 +57,8 @@ typedef struct s_linkcell { /* moldyn schedule structure */ typedef struct s_moldyn_schedule { - int content_count; + int count; + int total_sched; int *runs; double *tau; int (*hook)(void *moldyn,void *hook); diff --git a/run b/run index 3a15543..5a3a768 100755 --- a/run +++ b/run @@ -1,6 +1,6 @@ mkdir -p saves video ./clean $1 -./sic $1 +./sic $@ if [ "$?" == "0" ]; then #./perms if [ "$1" ] ; then diff --git a/sic.c b/sic.c index f4cdd7d..35931e0 100644 --- a/sic.c +++ b/sic.c @@ -14,9 +14,8 @@ int main(int argc,char **argv) { /* check argv */ - if(argc!=2) { - printf("[sic] error: arg1 (vis/log/save location) "); - printf("must be given!\n"); + if(argc!=3) { + printf("[sic] usage: %s \n",argv[0]); return -1; } @@ -32,7 +31,7 @@ int main(int argc,char **argv) { double tau; /* testing location & velocity vector */ - //t_3dvec r,v; + t_3dvec r,v; /* values */ tau=1.0e-15; /* delta t = 1 fs */ @@ -119,29 +118,30 @@ int main(int argc,char **argv) { moldyn_bc_check(&md); /* testing configuration */ - //r.x=2.45/2; v.x=0; + //r.x=2.95/2; v.x=0; //r.y=0; v.y=0; //r.z=0; v.z=0; //add_atom(&md,SI,M_SI,0, - // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP,//|ATOM_ATTR_HB, // &r,&v); - //r.x=-2.45/2; v.x=0; + //r.x=-2.95/2; v.x=0; //r.y=0; v.y=0; //r.z=0; v.z=0; //add_atom(&md,SI,M_SI,0, - // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP,//|ATOM_ATTR_HB, // &r,&v); /* setting a nearest neighbour distance for the moldyn checks */ set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */ /* set temperature */ - printf("[sic] setting temperature\n"); + printf("[sic] setting temperature -> %f\n",273+atof(argv[2])); //set_temperature(&md,273.0+1410.0); - set_temperature(&md,273.0+450.0); + //set_temperature(&md,273.0+450.0); //set_temperature(&md,273.0); //set_temperature(&md,1.0); //set_temperature(&md,0.0); + set_temperature(&md,atof(argv[2])+273.0); /* set pressure */ printf("[sic] setting pressure\n"); @@ -159,13 +159,17 @@ int main(int argc,char **argv) { /* create the simulation schedule */ printf("[sic] adding schedule\n"); - moldyn_add_schedule(&md,20001,1.0); + moldyn_add_schedule(&md,20000,.1); + moldyn_add_schedule(&md,10000,.2); + moldyn_add_schedule(&md,6667,.3); + moldyn_add_schedule(&md,5000,.4); + moldyn_add_schedule(&md,4001,.5); /* activate logging */ printf("[sic] activate logging\n"); moldyn_set_log_dir(&md,argv[1]); - moldyn_set_log(&md,LOG_TOTAL_ENERGY,200); - moldyn_set_log(&md,VISUAL_STEP,200); + moldyn_set_log(&md,LOG_TOTAL_ENERGY,100); + moldyn_set_log(&md,VISUAL_STEP,100); /* * let's do the actual md algorithm now -- 2.20.1