added msd calc + 11x constraints as ifdef

diff --git a/moldyn.c b/moldyn.c
index 15dcede..bd29afc 100644 (file)
--- 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;i<count;i++) {
                /* check whether fixed atom */
@@ -2288,6 +2318,30 @@ int velocity_verlet(t_moldyn *moldyn) {
                        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
@@ -2296,6 +2350,30 @@ int velocity_verlet(t_moldyn *moldyn) {
                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]));
@@ -2328,6 +2406,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;i<count;i++) {
                /* check whether fixed atom */
@@ -3197,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;

diff --git a/moldyn.h b/moldyn.h
index bdb3758..a1866a9 100644 (file)
--- 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 (file)
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 <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+
+#include "moldyn.h"
+
+int usage(char *prog) {
+
+       printf("\nusage:\n");
+       printf("  %s <save file>\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;
+}
+