forgot the lattice constant

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 7698c24..d3d63eb 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -78,6 +78,10 @@ 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
        return 0;
 }
 
@@ -258,7 +262,7 @@ int set_potential(t_moldyn *moldyn,u8 type) {
 
        switch(type) {
                case MOLDYN_POTENTIAL_TM:
-                       moldyn->func1b=tersoff_mult_1bp;
+                       //moldyn->func1b=tersoff_mult_1bp;
                        moldyn->func3b_j1=tersoff_mult_3bp_j1;
                        moldyn->func3b_k1=tersoff_mult_3bp_k1;
                        moldyn->func3b_j2=tersoff_mult_3bp_j2;
@@ -802,6 +806,8 @@ 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;\
@@ -2071,11 +2077,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
                temperature_calc(moldyn);
                virial_sum(moldyn);
                pressure_calc(moldyn);
-               /*
+#ifdef PDEBUG  
                thermodynamic_pressure_calc(moldyn);
                printf("\n\nDEBUG: numeric pressure calc: %f\n\n",
                       moldyn->tp/BAR);
-               */
+#endif
 
                /* calculate fluctuations + averages */
                average_and_fluctuation_calc(moldyn);
@@ -2182,15 +2188,20 @@ int moldyn_integrate(t_moldyn *moldyn) {
                }
 
                /* display progress */
+#ifndef PDEBUG
                if(!(i%10)) {
+#endif
                        /* get current time */
                        gettimeofday(&t2,NULL);
 
 printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
        sched->count,i,moldyn->total_steps,
        moldyn->t,moldyn->t_avg,
+#ifndef PDEBUG
        moldyn->p/BAR,moldyn->p_avg/BAR,
-       //moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
+#else
+       moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
+#endif
        moldyn->volume,
        (int)(t2.tv_sec-t1.tv_sec));
 
@@ -2198,7 +2209,9 @@ printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
 
                        /* copy over time */
                        t1=t2;
+#ifndef PDEBUG
                }
+#endif
 
                /* increase absolute time */
                moldyn->time+=moldyn->tau;
@@ -2219,6 +2232,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;
 }
 
@@ -2236,6 +2265,12 @@ int velocity_verlet(t_moldyn *moldyn) {
        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
        for(i=0;i<count;i++) {
                /* check whether fixed atom */
                if(atom[i].attr&ATOM_ATTR_FP)
@@ -2243,9 +2278,20 @@ 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
                v3_add(&(atom[i].r),&(atom[i].r),&delta);
                v3_scale(&delta,&(atom[i].f),h*tau_square);
+#ifdef CONSTRAINT_110_5832
+               if(i==5832) {
+                       delta.y=-delta.x;
+               }
+#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));
 
                /* velocities [actually v(t+tau/2)] */
@@ -2261,11 +2307,21 @@ int velocity_verlet(t_moldyn *moldyn) {
 
        /* forces depending on chosen potential */
 #ifndef ALBE_FAST
-       potential_force_calc(moldyn);
+       // if albe, use albe force calc routine
+       //if(moldyn->func3b_j1==albe_mult_3bp_j1)
+       //      albe_potential_force_calc(moldyn);
+       //else
+               potential_force_calc(moldyn);
 #else
        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
        for(i=0;i<count;i++) {
                /* check whether fixed atom */
                if(atom[i].attr&ATOM_ATTR_FP)
@@ -2713,6 +2769,51 @@ int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
        return 0;
 }
         
+int check_per_bound_and_care_for_pbc(t_moldyn *moldyn,t_atom *a) {
+       
+       double x,y,z;
+       t_3dvec *dim;
+
+       dim=&(moldyn->dim);
+
+       x=dim->x/2;
+       y=dim->y/2;
+       z=dim->z/2;
+
+       if(moldyn->status&MOLDYN_STAT_PBX) {
+               if(a->r.x>=x) {
+                       a->pbc[0]+=1;
+                       a->r.x-=dim->x;
+               }
+               else if(-a->r.x>x) {
+                       a->pbc[0]-=1;
+                       a->r.x+=dim->x;
+               }
+       }
+       if(moldyn->status&MOLDYN_STAT_PBY) {
+               if(a->r.y>=y) {
+                       a->pbc[1]+=1;
+                       a->r.y-=dim->y;
+               }
+               else if(-a->r.y>y) {
+                       a->pbc[1]-=1;
+                       a->r.y+=dim->y;
+               }
+       }
+       if(moldyn->status&MOLDYN_STAT_PBZ) {
+               if(a->r.z>=z) {
+                       a->pbc[2]+=1;
+                       a->r.z-=dim->z;
+               }
+               else if(-a->r.z>z) {
+                       a->pbc[2]-=1;
+                       a->r.z+=dim->z;
+               }
+       }
+
+       return 0;
+}
+        
 /*
  * debugging / critical check functions
  */
@@ -2803,7 +2904,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);
@@ -3048,6 +3149,7 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
        int i;
        t_atom *atom;
        t_3dvec dist;
+       t_3dvec final_r;
        double d2;
        int a_cnt;
        int b_cnt;
@@ -3061,8 +3163,12 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
 
        for(i=0;i<moldyn->count;i++) {
 
-               v3_sub(&dist,&(atom[i].r),&(atom[i].r_0));
-               check_per_bound(moldyn,&dist);
+               /* care for pb crossing */
+               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) {