completed crt implementation

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 15dcede..a6f5c45 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -55,6 +55,11 @@ pthread_mutex_t *amutex;
 pthread_mutex_t emutex;
 #endif
 
+/* fully constrained relaxation technique - global pointers */
+u8 crtt;
+u8 *constraints;
+double *trafo_angle;
+
 /*
  * the moldyn functions
  */
@@ -78,10 +83,9 @@ 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
+       if(crtt)
+               printf("USING CRT\n");
+
        return 0;
 }
 
@@ -2256,6 +2260,43 @@ printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
        return 0;
 }
 
+/* basis trafo */
+
+#define FORWARD                0
+#define BACKWARD       1
+
+int basis_trafo(t_3dvec *r,u8 dir,double z,double x) {
+
+       t_3dvec tmp;
+
+       if(dir==FORWARD) {
+               if(z!=0.0) {
+                       v3_copy(&tmp,r);
+                       r->x=cos(z)*tmp.x-sin(z)*tmp.y;
+                       r->y=sin(z)*tmp.x+cos(z)*tmp.y;
+               }
+               if(x!=0.0) {
+                       v3_copy(&tmp,r);
+                       r->y=cos(x)*tmp.y-sin(x)*tmp.z;
+                       r->z=sin(x)*tmp.y+cos(x)*tmp.z;
+               }
+       }
+       else {
+               if(x!=0.0) {
+                       v3_copy(&tmp,r);
+                       r->y=cos(-x)*tmp.y-sin(-x)*tmp.z;
+                       r->z=sin(-x)*tmp.y+cos(-x)*tmp.z;
+               }
+               if(z!=0.0) {
+                       v3_copy(&tmp,r);
+                       r->x=cos(-z)*tmp.x-sin(-z)*tmp.y;
+                       r->y=sin(-z)*tmp.x+cos(-z)*tmp.y;
+               }
+       }
+
+       return 0;
+}
+
 /* velocity verlet */
 
 int velocity_verlet(t_moldyn *moldyn) {
@@ -2270,33 +2311,34 @@ 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)
                        continue;
+
                /* 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;
+
+               /* constraint relaxation */
+               if(crtt) {
+                       basis_trafo(&(atom[i].f),FORWARD,
+                                   trafo_angle[2*i],trafo_angle[2*i+1]);
+                       if(constraints[3*i])
+                               atom[i].f.x=0;
+                       if(constraints[3*i+1])
+                               atom[i].f.y=0;
+                       if(constraints[3*i+2])
+                               atom[i].f.z=0;
+                       basis_trafo(&(atom[i].f),BACKWARD,
+                                   trafo_angle[2*i],trafo_angle[2*i+1]);
                }
-#endif
+
 #ifndef QUENCH
+               v3_scale(&delta,&(atom[i].v),tau);
                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
                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));
@@ -2323,16 +2365,25 @@ 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
        for(i=0;i<count;i++) {
                /* check whether fixed atom */
                if(atom[i].attr&ATOM_ATTR_FP)
                        continue;
+
+               /* constraint relaxation */
+               if(crtt) {
+                       basis_trafo(&(atom[i].f),FORWARD,
+                                   trafo_angle[2*i],trafo_angle[2*i+1]);
+                       if(constraints[3*i])
+                               atom[i].f.x=0;
+                       if(constraints[3*i+1])
+                               atom[i].f.y=0;
+                       if(constraints[3*i+2])
+                               atom[i].f.z=0;
+                       basis_trafo(&(atom[i].f),BACKWARD,
+                                   trafo_angle[2*i],trafo_angle[2*i+1]);
+               }
+
                /* again velocities [actually v(t+tau)] */
                v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
                v3_add(&(atom[i].v),&(atom[i].v),&delta);
@@ -3197,6 +3248,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;