X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=415581c475fef1e8bb85b4fb8a1595135de3a88f;hb=0d2926a23ea7bbf7654da88770c827ad041357a5;hp=2fdf4c2e0b8a7bf47a3f3fb7dc76d235465c48e6;hpb=5b64e8a2762ea3248792a3c8cb57add0ecb5442c;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 2fdf4c2..415581c 100644 --- a/moldyn.c +++ b/moldyn.c @@ -13,6 +13,8 @@ #include #include #include +#include +#include #include #include "moldyn.h" @@ -24,6 +26,9 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) { memset(moldyn,0,sizeof(t_moldyn)); + moldyn->argc=argc; + moldyn->args=argv; + rand_init(&(moldyn->random),NULL,1); moldyn->random.status|=RAND_STAT_VERBOSE; @@ -510,6 +515,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, atom[ret].brand=brand; atom[ret].tag=count+ret; check_per_bound(moldyn,&(atom[ret].r)); + atom[ret].r_0=atom[ret].r; } /* update total system mass */ @@ -518,6 +524,40 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, return ret; } +int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, + t_3dvec *r,t_3dvec *v) { + + t_atom *atom; + void *ptr; + int count; + + atom=moldyn->atom; + count=(moldyn->count)++; + + ptr=realloc(atom,(count+1)*sizeof(t_atom)); + if(!ptr) { + perror("[moldyn] realloc (add atom)"); + return -1; + } + moldyn->atom=ptr; + + atom=moldyn->atom; + atom[count].r=*r; + atom[count].v=*v; + atom[count].element=element; + atom[count].mass=mass; + atom[count].brand=brand; + atom[count].tag=count; + atom[count].attr=attr; + check_per_bound(moldyn,&(atom[count].r)); + atom[count].r_0=atom[count].r; + + /* update total system mass */ + total_mass_calc(moldyn); + + return 0; +} + /* cubic init */ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { @@ -630,38 +670,6 @@ int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { return count; } -int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, - t_3dvec *r,t_3dvec *v) { - - t_atom *atom; - void *ptr; - int count; - - atom=moldyn->atom; - count=(moldyn->count)++; - - ptr=realloc(atom,(count+1)*sizeof(t_atom)); - if(!ptr) { - perror("[moldyn] realloc (add atom)"); - return -1; - } - moldyn->atom=ptr; - - atom=moldyn->atom; - atom[count].r=*r; - atom[count].v=*v; - atom[count].element=element; - atom[count].mass=mass; - atom[count].brand=brand; - atom[count].tag=count; - atom[count].attr=attr; - - /* update total system mass */ - total_mass_calc(moldyn); - - return 0; -} - int destroy_atoms(t_moldyn *moldyn) { if(moldyn->atom) free(moldyn->atom); @@ -1097,8 +1105,10 @@ double e_kin_calc(t_moldyn *moldyn) { atom=moldyn->atom; moldyn->ekin=0.0; - for(i=0;icount;i++) - moldyn->ekin+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v)); + for(i=0;icount;i++) { + atom[i].ekin=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v)); + moldyn->ekin+=atom[i].ekin; + } return moldyn->ekin; } @@ -1340,6 +1350,7 @@ int moldyn_integrate(t_moldyn *moldyn) { char dir[128]; double ds; double energy_scale; + struct timeval t1,t2; //double tp; sched=&(moldyn->schedule); @@ -1360,8 +1371,8 @@ int moldyn_integrate(t_moldyn *moldyn) { moldyn->tau_square=moldyn->tau*moldyn->tau; moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff; - /* energy scaling factor */ - energy_scale=moldyn->count*EV; + /* get current time */ + gettimeofday(&t1,NULL); /* calculate initial forces */ potential_force_calc(moldyn); @@ -1399,6 +1410,9 @@ return 0; moldyn->tau_square=moldyn->tau*moldyn->tau; moldyn->time_steps=sched->runs[sched->count]; + /* energy scaling factor (might change!) */ + energy_scale=moldyn->count*EV; + /* integration according to schedule */ for(i=0;itime_steps;i++) { @@ -1421,7 +1435,7 @@ return 0; /* check for log & visualization */ if(e) { - if(!(i%e)) + if(!(moldyn->total_steps%e)) dprintf(moldyn->efd, "%f %f %f %f\n", moldyn->time,moldyn->ekin/energy_scale, @@ -1429,7 +1443,7 @@ return 0; get_total_energy(moldyn)/energy_scale); } if(m) { - if(!(i%m)) { + if(!(moldyn->total_steps%m)) { momentum=get_total_p(moldyn); dprintf(moldyn->mfd, "%f %f %f %f %f\n",moldyn->time, @@ -1438,7 +1452,7 @@ return 0; } } if(p) { - if(!(i%p)) { + if(!(moldyn->total_steps%p)) { dprintf(moldyn->pfd, "%f %f %f %f %f\n",moldyn->time, moldyn->p/BAR,moldyn->p_avg/BAR, @@ -1446,17 +1460,18 @@ return 0; } } if(t) { - if(!(i%t)) { + if(!(moldyn->total_steps%t)) { dprintf(moldyn->tfd, "%f %f %f\n", moldyn->time,moldyn->t,moldyn->t_avg); } } if(s) { - if(!(i%s)) { + if(!(moldyn->total_steps%s)) { snprintf(dir,128,"%s/s-%07.f.save", moldyn->vlsdir,moldyn->time); - fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT); + fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT, + S_IRUSR|S_IWUSR); if(fd<0) perror("[moldyn] save fd open"); else { write(fd,moldyn,sizeof(t_moldyn)); @@ -1467,20 +1482,27 @@ return 0; } } if(v) { - if(!(i%v)) { + if(!(moldyn->total_steps%v)) { visual_atoms(&(moldyn->vis),moldyn->time, moldyn->atom,moldyn->count); } } /* display progress */ - if(!(i%10)) { - printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f", + if(!(moldyn->total_steps%10)) { + /* get current time */ + gettimeofday(&t2,NULL); + + printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)", sched->count,i, moldyn->t,moldyn->t_avg, - moldyn->p_avg/BAR,moldyn->p/BAR, - moldyn->volume); + moldyn->p_avg/BAR,moldyn->gp_avg/BAR, + moldyn->volume, + (int)(t2.tv_sec-t1.tv_sec)); fflush(stdout); + + /* copy over time */ + t1=t2; } /* increase absolute time */ @@ -1490,13 +1512,11 @@ return 0; } /* check for hooks */ - if(sched->count+1total_sched) { - if(sched->hook) { - printf("\n ## schedule hook %d/%d start ##\n", - sched->count+1,sched->total_sched-1); - sched->hook(moldyn,sched->hook_params); - printf(" ## schedule hook end ##\n"); - } + if(sched->hook) { + printf("\n ## schedule hook %d/%d start ##\n", + sched->count+1,sched->total_sched-1); + sched->hook(moldyn,sched->hook_params); + printf(" ## schedule hook end ##\n"); } /* increase the schedule counter */ @@ -1917,6 +1937,17 @@ int moldyn_bc_check(t_moldyn *moldyn) { return 0; } +/* + * restore function + */ + +int moldyn_load(t_moldyn *moldyn) { + + // later ... + + return 0; +} + /* * post processing functions */ @@ -1939,3 +1970,11 @@ int get_line(int fd,char *line,int max) { } } +int analyze_bonds(t_moldyn *moldyn) { + + + + + return 0; +} +