X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=3c4d87399dcaceffb42f6d91ffe97a58bd8e841f;hb=c40d54eb3e319b17b2f6174c4eddcfd6ee3a407b;hp=cd2a80309c6cf9de56e4b23132debbe9e46b02cd;hpb=0fbd532d8ddce9848df592ed586ffba439369284;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index cd2a803..3c4d873 100644 --- a/moldyn.c +++ b/moldyn.c @@ -13,17 +13,74 @@ #include #include #include +#include +#include #include #include "moldyn.h" #include "report/report.h" +/* + * global variables, pse and atom colors (only needed here) + */ + +static char *pse_name[]={ + "*", + "H", + "He", + "Li", + "Be", + "B", + "C", + "N", + "O", + "F", + "Ne", + "Na", + "Mg", + "Al", + "Si", + "P", + "S", + "Cl", + "Ar", +}; + +static char *pse_col[]={ + "*", + "White", + "He", + "Li", + "Be", + "B", + "Gray", + "N", + "Blue", + "F", + "Ne", + "Na", + "Mg", + "Al", + "Yellow", + "P", + "S", + "Cl", + "Ar", +}; + +/* + * the moldyn functions + */ + int moldyn_init(t_moldyn *moldyn,int argc,char **argv) { printf("[moldyn] init\n"); memset(moldyn,0,sizeof(t_moldyn)); + moldyn->argc=argc; + moldyn->args=argv; + rand_init(&(moldyn->random),NULL,1); moldyn->random.status|=RAND_STAT_VERBOSE; @@ -218,6 +275,14 @@ int set_potential_params(t_moldyn *moldyn,void *params) { return 0; } +int set_avg_skip(t_moldyn *moldyn,int skip) { + + printf("[moldyn] skip %d steps before starting average calc\n",skip); + moldyn->avg_skip=skip; + + return 0; +} + int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) { strncpy(moldyn->vlsdir,dir,127); @@ -414,7 +479,6 @@ int moldyn_log_shutdown(t_moldyn *moldyn) { moldyn->vlsdir); system(sc); } - if(&(moldyn->vis)) visual_tini(&(moldyn->vis)); return 0; } @@ -502,6 +566,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 */ @@ -510,6 +575,69 @@ 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; +} + +int del_atom(t_moldyn *moldyn,int tag) { + + t_atom *new,*old; + int cnt; + + old=moldyn->atom; + + new=(t_atom *)malloc((moldyn->count-1)*sizeof(t_atom)); + if(!new) { + perror("[moldyn]malloc (del atom)"); + return -1; + } + + for(cnt=0;cntcount;cnt++) { + new[cnt-1]=old[cnt]; + new[cnt-1].tag=cnt-1; + } + + moldyn->count-=1; + moldyn->atom=new; + + free(old); + + return 0; +} + /* cubic init */ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { @@ -622,38 +750,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); @@ -728,8 +824,6 @@ double temperature_calc(t_moldyn *moldyn) { /* assume up to date kinetic energy, which is 3/2 N k_B T */ moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count); - moldyn->t_sum+=moldyn->t; - moldyn->mean_t=moldyn->t_sum/moldyn->total_steps; return moldyn->t; } @@ -804,64 +898,88 @@ double ideal_gas_law_pressure(t_moldyn *moldyn) { return p; } -double pressure_calc(t_moldyn *moldyn) { +double virial_sum(t_moldyn *moldyn) { int i; double v; t_virial *virial; + /* virial (sum over atom virials) */ + v=0.0; + for(i=0;icount;i++) { + virial=&(moldyn->atom[i].virial); + v+=(virial->xx+virial->yy+virial->zz); + } + moldyn->virial=v; + + /* global virial (absolute coordinates) */ + virial=&(moldyn->gvir); + moldyn->gv=virial->xx+virial->yy+virial->zz; + + return moldyn->virial; +} + +double pressure_calc(t_moldyn *moldyn) { + /* * PV = NkT + - * W = 1/3 sum_i f_i r_i + * with W = 1/3 sum_i f_i r_i (- skipped!) * virial = sum_i f_i r_i * * => P = (2 Ekin + virial) / (3V) */ - v=0.0; - for(i=0;icount;i++) { - virial=&(moldyn->atom[i].virial); - v+=(virial->xx+virial->yy+virial->zz); - } - - /* virial sum and mean virial */ - moldyn->virial_sum+=v; - moldyn->mean_v=moldyn->virial_sum/moldyn->total_steps; + /* assume up to date virial & up to date kinetic energy */ - /* assume up to date kinetic energy */ - moldyn->p=2.0*moldyn->ekin+moldyn->mean_v; + /* pressure (atom virials) */ + moldyn->p=2.0*moldyn->ekin+moldyn->virial; moldyn->p/=(3.0*moldyn->volume); - moldyn->p_sum+=moldyn->p; - moldyn->mean_p=moldyn->p_sum/moldyn->total_steps; - /* pressure from 'absolute coordinates' virial */ - virial=&(moldyn->virial); - v=virial->xx+virial->yy+virial->zz; - moldyn->gp=2.0*moldyn->ekin+v; + /* pressure (absolute coordinates) */ + moldyn->gp=2.0*moldyn->ekin+moldyn->gv; moldyn->gp/=(3.0*moldyn->volume); - moldyn->gp_sum+=moldyn->gp; - moldyn->mean_gp=moldyn->gp_sum/moldyn->total_steps; return moldyn->p; } -int energy_fluctuation_calc(t_moldyn *moldyn) { +int average_and_fluctuation_calc(t_moldyn *moldyn) { + + if(moldyn->total_stepsavg_skip) + return 0; - /* assume up to date energies */ + int denom=moldyn->total_steps+1-moldyn->avg_skip; + + /* assume up to date energies, temperature, pressure etc */ /* kinetic energy */ moldyn->k_sum+=moldyn->ekin; moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin); - moldyn->k_mean=moldyn->k_sum/moldyn->total_steps; - moldyn->k2_mean=moldyn->k2_sum/moldyn->total_steps; - moldyn->dk2_mean=moldyn->k2_mean-(moldyn->k_mean*moldyn->k_mean); + moldyn->k_avg=moldyn->k_sum/denom; + moldyn->k2_avg=moldyn->k2_sum/denom; + moldyn->dk2_avg=moldyn->k2_avg-(moldyn->k_avg*moldyn->k_avg); /* potential energy */ moldyn->v_sum+=moldyn->energy; moldyn->v2_sum+=(moldyn->energy*moldyn->energy); - moldyn->v_mean=moldyn->v_sum/moldyn->total_steps; - moldyn->v2_mean=moldyn->v2_sum/moldyn->total_steps; - moldyn->dv2_mean=moldyn->v2_mean-(moldyn->v_mean*moldyn->v_mean); + moldyn->v_avg=moldyn->v_sum/denom; + moldyn->v2_avg=moldyn->v2_sum/denom; + moldyn->dv2_avg=moldyn->v2_avg-(moldyn->v_avg*moldyn->v_avg); + + /* temperature */ + moldyn->t_sum+=moldyn->t; + moldyn->t_avg=moldyn->t_sum/denom; + + /* virial */ + moldyn->virial_sum+=moldyn->virial; + moldyn->virial_avg=moldyn->virial_sum/denom; + moldyn->gv_sum+=moldyn->gv; + moldyn->gv_avg=moldyn->gv_sum/denom; + + /* pressure */ + moldyn->p_sum+=moldyn->p; + moldyn->p_avg=moldyn->p_sum/denom; + moldyn->gp_sum+=moldyn->gp; + moldyn->gp_avg=moldyn->gp_sum/denom; return 0; } @@ -870,10 +988,14 @@ int get_heat_capacity(t_moldyn *moldyn) { double temp2,ighc; + /* averages needed for heat capacity calc */ + if(moldyn->total_stepsavg_skip) + return 0; + /* (temperature average)^2 */ - temp2=moldyn->mean_t*moldyn->mean_t; + temp2=moldyn->t_avg*moldyn->t_avg; printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n", - moldyn->mean_t); + moldyn->t_avg); /* ideal gas contribution */ ighc=3.0*moldyn->count*K_BOLTZMANN/2.0; @@ -881,15 +1003,16 @@ int get_heat_capacity(t_moldyn *moldyn) { ighc/moldyn->mass*KILOGRAM/JOULE); /* specific heat for nvt ensemble */ - moldyn->c_v_nvt=moldyn->dv2_mean/(K_BOLTZMANN*temp2)+ighc; + moldyn->c_v_nvt=moldyn->dv2_avg/(K_BOLTZMANN*temp2)+ighc; moldyn->c_v_nvt/=moldyn->mass; /* specific heat for nve ensemble */ - moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_mean/(ighc*K_BOLTZMANN*temp2))); + moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_avg/(ighc*K_BOLTZMANN*temp2))); moldyn->c_v_nve/=moldyn->mass; printf(" NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE); printf(" NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE); +printf(" --> sim: %f experimental: %f\n",moldyn->dv2_avg,1.5*moldyn->count*K_B2*moldyn->t_avg*moldyn->t_avg*(1.0-1.5*moldyn->count*K_BOLTZMANN/(700*moldyn->mass*JOULE/KILOGRAM))); return 0; } @@ -1062,8 +1185,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; } @@ -1305,6 +1430,7 @@ int moldyn_integrate(t_moldyn *moldyn) { char dir[128]; double ds; double energy_scale; + struct timeval t1,t2; //double tp; sched=&(moldyn->schedule); @@ -1325,8 +1451,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); @@ -1356,13 +1482,17 @@ return 0; printf("[moldyn] integration start, go get a coffee ...\n"); /* executing the schedule */ - for(sched->count=0;sched->counttotal_sched;sched->count++) { + sched->count=0; + while(sched->counttotal_sched) { /* setting amount of runs and finite time step size */ moldyn->tau=sched->tau[sched->count]; 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++) { @@ -1373,10 +1503,9 @@ return 0; /* calculate kinetic energy, temperature and pressure */ e_kin_calc(moldyn); temperature_calc(moldyn); + virial_sum(moldyn); pressure_calc(moldyn); - energy_fluctuation_calc(moldyn); - //tp=thermodynamic_pressure_calc(moldyn); -//printf("thermodynamic p: %f\n",thermodynamic_pressure_calc(moldyn)/BAR); + average_and_fluctuation_calc(moldyn); /* p/t scaling */ if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT)) @@ -1386,7 +1515,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, @@ -1394,7 +1523,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, @@ -1403,25 +1532,26 @@ 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->mean_p/BAR, - moldyn->gp/BAR,moldyn->mean_gp/BAR); + moldyn->p/BAR,moldyn->p_avg/BAR, + moldyn->gp/BAR,moldyn->gp_avg/BAR); } } if(t) { - if(!(i%t)) { + if(!(moldyn->total_steps%t)) { dprintf(moldyn->tfd, "%f %f %f\n", - moldyn->time,moldyn->t,moldyn->mean_t); + 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)); @@ -1432,23 +1562,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: %f, P: %f %f V: %f", - sched->count,i, - moldyn->mean_t, - moldyn->mean_p/BAR, - moldyn->mean_gp/BAR, - moldyn->volume); - fflush(stdout); -printf("\n"); -get_heat_capacity(moldyn); + 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->gp_avg/BAR, + moldyn->volume, + (int)(t2.tv_sec-t1.tv_sec)); + fflush(stdout); + + /* copy over time */ + t1=t2; } /* increase absolute time */ @@ -1458,12 +1592,15 @@ get_heat_capacity(moldyn); } /* check for hooks */ - if(sched->count+1total_sched) - if(sched->hook) - sched->hook(moldyn,sched->hook_params); + 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"); + } - /* get a new info line */ - printf("\n"); + /* increase the schedule counter */ + sched->count+=1; } @@ -1542,7 +1679,7 @@ int potential_force_calc(t_moldyn *moldyn) { moldyn->energy=0.0; /* reset global virial */ - memset(&(moldyn->virial),0,sizeof(t_virial)); + memset(&(moldyn->gvir),0,sizeof(t_virial)); /* reset force, site energy and virial of every atom */ for(i=0;ivirial.xx+=moldyn->atom[i].r.x*moldyn->atom[i].f.x; - moldyn->virial.yy+=moldyn->atom[i].r.y*moldyn->atom[i].f.y; - moldyn->virial.zz+=moldyn->atom[i].r.z*moldyn->atom[i].f.z; - moldyn->virial.xy+=moldyn->atom[i].r.y*moldyn->atom[i].f.x; - moldyn->virial.xz+=moldyn->atom[i].r.z*moldyn->atom[i].f.x; - moldyn->virial.yz+=moldyn->atom[i].r.z*moldyn->atom[i].f.y; + moldyn->gvir.xx+=moldyn->atom[i].r.x*moldyn->atom[i].f.x; + moldyn->gvir.yy+=moldyn->atom[i].r.y*moldyn->atom[i].f.y; + moldyn->gvir.zz+=moldyn->atom[i].r.z*moldyn->atom[i].f.z; + moldyn->gvir.xy+=moldyn->atom[i].r.y*moldyn->atom[i].f.x; + moldyn->gvir.xz+=moldyn->atom[i].r.z*moldyn->atom[i].f.x; + moldyn->gvir.yz+=moldyn->atom[i].r.z*moldyn->atom[i].f.y; } return 0; @@ -1880,6 +2017,17 @@ int moldyn_bc_check(t_moldyn *moldyn) { return 0; } +/* + * restore function + */ + +int moldyn_load(t_moldyn *moldyn) { + + // later ... + + return 0; +} + /* * post processing functions */ @@ -1902,3 +2050,101 @@ int get_line(int fd,char *line,int max) { } } +int analyze_bonds(t_moldyn *moldyn) { + + + + + return 0; +} + +/* + * visualization code + */ + +int visual_init(t_visual *v,char *filebase) { + + char file[128+8]; + + strncpy(v->fb,filebase,128); + memset(file,0,128+8); + + return 0; +} + +int visual_atoms(t_visual *v,double time,t_atom *atom,int n) { + + int i,fd; + char file[128+64]; + t_3dvec dim; + double help; + + dim.x=v->dim.x; + dim.y=v->dim.y; + dim.z=v->dim.z; + + help=(dim.x+dim.y); + + sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,time); + fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR); + if(fd<0) { + perror("open visual save file fd"); + return -1; + } + + /* write the actual data file */ + dprintf(fd,"# [P] %d %07.f <%f,%f,%f>\n", + n,time,help/40.0,help/40.0,-0.8*help); + for(i=0;i