From: hackbard Date: Fri, 11 Sep 2009 16:31:39 +0000 (+0200) Subject: added msd calc + 11x constraints as ifdef X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=cba0d347201db1c99dc23736658ed78ed12ad5bd;p=physik%2Fposic.git added msd calc + 11x constraints as ifdef --- diff --git a/moldyn.c b/moldyn.c index 15dcede..bd29afc 100644 --- a/moldyn.c +++ b/moldyn.c @@ -81,6 +81,10 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) { #ifdef CONSTRAINT_110_5832 printf("\n\n\nWARNING! WARNING! WARNING!\n\n\n"); printf("\n\n\n!! -- constraints enabled -- !!\n\n\n"); +#endif +#ifdef CONSTRAINT_11X_5832 + printf("\n\n\nWARNING! WARNING! WARNING!\n\n\n"); + printf("\n\n\n!! -- constraints enabled -- !!\n\n\n"); #endif return 0; } @@ -2264,6 +2268,10 @@ int velocity_verlet(t_moldyn *moldyn) { double tau,tau_square,h; t_3dvec delta; t_atom *atom; +#ifdef CONSTRAINT_11X_5832 + double xt,yt,zt; + double xtt,ytt,ztt; +#endif atom=moldyn->atom; count=moldyn->count; @@ -2275,6 +2283,28 @@ int velocity_verlet(t_moldyn *moldyn) { atom[5832].f.x=0.5*(atom[5832].f.x-atom[5832].f.y); atom[5832].f.y=-atom[5832].f.x; } +#endif +#ifdef CONSTRAINT_11X_5832 + if(count==5833) { + // second trafo + xt=atom[5832].f.x; + yt=atom[5832].f.y*cos(-0.16935129)-atom[5832].f.z*sin(-0.16935129); + zt=atom[5832].f.y*sin(-0.16935129)+atom[5832].f.z*cos(-0.16935129); + // first trafo + xtt=xt*cos(M_PI/4.0)-yt*sin(M_PI/4.0); + ytt=xt*sin(M_PI/4.0)+yt*sin(M_PI/4.0); + ztt=zt; + // apply constraints + ytt=0.0; + // first trafo backwards + xt=xtt*cos(M_PI/4.0)+ytt*sin(M_PI/4.0); + yt=-xtt*sin(M_PI/4.0)+ytt*sin(M_PI/4.0); + zt=ztt; + // second trafo backwards + atom[5832].f.x=xt; + atom[5832].f.y=yt*cos(-0.16935129)+zt*sin(-0.16935129); + atom[5832].f.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129); + } #endif for(i=0;iatom; + msd[0]=0; + msd[1]=0; + msd[2]=0; + a_cnt=0; + b_cnt=0; + + for(i=0;icount;i++) { + + /* care for pb crossing */ + if(atom[i].pbc[0]|atom[i].pbc[1]|atom[i].pbc[2]) { + printf("[moldyn] msd pb crossings for atom %d\n",i); + printf(" x: %d y: %d z: %d\n", + atom[i].pbc[0],atom[i].pbc[1],atom[i].pbc[2]); + } + final_r.x=atom[i].r.x+atom[i].pbc[0]*moldyn->dim.x; + final_r.y=atom[i].r.y+atom[i].pbc[1]*moldyn->dim.y; + final_r.z=atom[i].r.z+atom[i].pbc[2]*moldyn->dim.z; + /* calculate distance */ + v3_sub(&dist,&final_r,&(atom[i].r_0)); + d2=v3_absolute_square(&dist); + + if(atom[i].brand) { + b_cnt+=1; + msd[1]+=d2; + } + else { + a_cnt+=1; + msd[0]+=d2; + } + + msd[2]+=d2; + } + + msd[0]/=a_cnt; + msd[1]/=b_cnt; + msd[2]/=moldyn->count; + + return 0; +} + int bonding_analyze(t_moldyn *moldyn,double *cnt) { return 0; diff --git a/moldyn.h b/moldyn.h index bdb3758..a1866a9 100644 --- a/moldyn.h +++ b/moldyn.h @@ -523,6 +523,7 @@ int get_line(int fd,char *line,int max); int pair_correlation_init(t_moldyn *moldyn,double dr); int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc); +int calculate_msd(t_moldyn *moldyn,double *msd); int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom, t_atom *jtom,void *data,u8 bc); int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr); diff --git a/msd_calc.c b/msd_calc.c new file mode 100644 index 0000000..80388d8 --- /dev/null +++ b/msd_calc.c @@ -0,0 +1,54 @@ +/* + * calculation of mean square displacement + * + * author: frank.zirkelbach@physik.uni-augsburg.de + * + */ + +#define _GNU_SOURCE +#include +#include +#include +#include +#include +#include +#include + +#include "moldyn.h" + +int usage(char *prog) { + + printf("\nusage:\n"); + printf(" %s \n\n",prog); + + return -1; +} + +int main(int argc,char **argv) { + + t_moldyn moldyn; + int ret; + double msd[3]; + + if(argc!=2) { + usage(argv[0]); + return -1; + } + + memset(&moldyn,0,sizeof(t_moldyn)); + + printf("[msd calc] reading save file ...\n"); + ret=moldyn_read_save_file(&moldyn,argv[1]); + if(ret) { + printf("[msd calc] exit!\n"); + return ret; + } + + calculate_msd(&moldyn,msd); + + printf("MSD - %f ps: %.10f %.10f %.10f\n",moldyn.time, + msd[0],msd[1],msd[2]); + + return 0; +} +