some checkins

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 12cb49b..bd29afc 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -78,6 +78,14 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
        pthread_mutex_init(&emutex,NULL);
 #endif
 
+#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;
 }
 
@@ -516,7 +524,8 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
 
 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
                    u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin,
-                   t_part_params *p_params,t_defect_params *d_params) {
+                   t_part_params *p_params,t_defect_params *d_params,
+                   t_offset_params *o_params) {
 
        int new,count;
        int ret;
@@ -588,6 +597,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
 
        switch(type) {
                case CUBIC:
+                       if(o_params->use)
+                               v3_add(&orig,&orig,&(o_params->o));
                        set_nn_dist(moldyn,lc);
                        ret=cubic_init(a,b,c,lc,atom,&orig,p_params,d_params);
                        strcpy(name,"cubic");
@@ -595,6 +606,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
                case FCC:
                        if(!origin)
                                v3_scale(&orig,&orig,0.5);
+                       if(o_params->use)
+                               v3_add(&orig,&orig,&(o_params->o));
                        set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
                        ret=fcc_init(a,b,c,lc,atom,&orig,p_params,d_params);
                        strcpy(name,"fcc");
@@ -602,6 +615,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
                case DIAMOND:
                        if(!origin)
                                v3_scale(&orig,&orig,0.25);
+                       if(o_params->use)
+                               v3_add(&orig,&orig,&(o_params->o));
                        set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
                        ret=diamond_init(a,b,c,lc,atom,&orig,p_params,d_params);
                        strcpy(name,"diamond");
@@ -802,8 +817,6 @@ int del_atom(t_moldyn *moldyn,int tag) {
                        case DEFECT_STYPE_DB_Z:\
                                d_o.z=d_params->od;\
                                d_d.z=d_params->dd;\
-d_d.x=0.9;\
-d_d.y=0.9;\
                                break;\
                        case DEFECT_STYPE_DB_R:\
                                break;\
@@ -2228,6 +2241,22 @@ printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
 
        }
 
+       /* writing a final save file! */
+       if(s) {
+               snprintf(dir,128,"%s/s-final.save",moldyn->vlsdir);
+               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));
+                       write(fd,moldyn->atom,
+                             moldyn->count*sizeof(t_atom));
+               }
+               close(fd);
+       }
+       /* writing a final visual file! */
+       if(a)
+               visual_atoms(moldyn);
+
        return 0;
 }
 
@@ -2239,12 +2268,44 @@ 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;
        tau=moldyn->tau;
        tau_square=moldyn->tau_square;
 
+#ifdef CONSTRAINT_110_5832
+       if(count==5833) {
+               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;i<count;i++) {
                /* check whether fixed atom */
                if(atom[i].attr&ATOM_ATTR_FP)
@@ -2252,8 +2313,68 @@ int velocity_verlet(t_moldyn *moldyn) {
                /* new positions */
                h=0.5/atom[i].mass;
                v3_scale(&delta,&(atom[i].v),tau);
+#ifdef CONSTRAINT_110_5832
+               if(i==5832) {
+                       delta.y=-delta.x;
+               }
+#endif
+#ifdef JHDVHJDV
+#ifdef CONSTRAINT_11X_5832
+       if(i==5833) {
+               // second trafo
+               xt=delta.x;
+               yt=delta.y*cos(-0.16935129)-delta.z*sin(-0.16935129);
+               zt=delta.y*sin(-0.16935129)+delta.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
+               delta.x=xt;
+               delta.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
+               delta.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
+       }
+#endif
+#endif
+#ifndef QUENCH
                v3_add(&(atom[i].r),&(atom[i].r),&delta);
+#endif
                v3_scale(&delta,&(atom[i].f),h*tau_square);
+#ifdef CONSTRAINT_110_5832
+               if(i==5832) {
+                       delta.y=-delta.x;
+               }
+#endif
+#ifdef GHDJHBSJBSD
+#ifdef CONSTRAINT_11X_5832
+       if(i==5833) {
+               // second trafo
+               xt=delta.x;
+               yt=delta.y*cos(-0.16935129)-delta.z*sin(-0.16935129);
+               zt=delta.y*sin(-0.16935129)+delta.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
+               delta.x=xt;
+               delta.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
+               delta.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
+       }
+#endif
+#endif
                v3_add(&(atom[i].r),&(atom[i].r),&delta);
                //check_per_bound_and_care_for_pbc(moldyn,&(atom[i]));
                check_per_bound(moldyn,&(atom[i].r));
@@ -2280,6 +2401,34 @@ int velocity_verlet(t_moldyn *moldyn) {
        albe_potential_force_calc(moldyn);
 #endif
 
+#ifdef CONSTRAINT_110_5832
+       if(count==5833) {
+               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;i<count;i++) {
                /* check whether fixed atom */
                if(atom[i].attr&ATOM_ATTR_FP)
@@ -2862,7 +3011,7 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
        if(fsize!=sizeof(t_moldyn)+size) {
                corr=fsize-sizeof(t_moldyn)-size;
                printf("[moldyn] WARNING: lsf (illegal file size)\n");
-               printf("  moifying offset:\n");
+               printf("  modifying offset:\n");
                printf("  - current pos: %d\n",sizeof(t_moldyn));
                printf("  - atom size: %d\n",size);
                printf("  - file size: %d\n",fsize);
@@ -3148,6 +3297,57 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
        return 0;
 }
 
+int calculate_msd(t_moldyn *moldyn,double *msd) {
+
+       int i;
+       t_atom *atom;
+       t_3dvec dist;
+       t_3dvec final_r;
+       double d2;
+       int a_cnt;
+       int b_cnt;
+
+       atom=moldyn->atom;
+       msd[0]=0;
+       msd[1]=0;
+       msd[2]=0;
+       a_cnt=0;
+       b_cnt=0;
+
+       for(i=0;i<moldyn->count;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;