insert mixed atoms feature added (will be renamed an cleand up later)

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 988b467..eb3c4b0 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -30,99 +30,6 @@
 #include "potentials/tersoff.h"
 #endif
 
-
-/*
- * global variables, pse and atom colors (only needed here)
- */ 
-
-static char *pse_name[]={
-       "*",
-       "H",
-       "He",
-       "Li",
-       "Be",
-       "B",
-       "C",
-       "N",
-       "O",
-       "F",
-       "Ne",
-       "Na",
-       "Mg",
-       "Al",
-       "Si",
-       "P",
-       "S",
-       "Cl",
-       "Ar",
-};
-
-static char *pse_col[]={
-       "*",
-       "White",
-       "He",
-       "Li",
-       "Be",
-       "B",
-       "Gray",
-       "N",
-       "Blue",
-       "F",
-       "Ne",
-       "Na",
-       "Mg",
-       "Al",
-       "Yellow",
-       "P",
-       "S",
-       "Cl",
-       "Ar",
-};
-
-static double *pse_mass[]={
-       0,
-       0,
-       0,
-       0,
-       0,
-       0,
-       M_C,
-       0,
-       0,
-       0,
-       0,
-       0,
-       0,
-       0,
-       M_SI,
-       0,
-       0,
-       0,
-       0,
-};
-
-static double *pse_lc[]={
-       0,
-       0,
-       0,
-       0,
-       0,
-       0,
-       LC_C,
-       0,
-       0,
-       0,
-       0,
-       0,
-       0,
-       0,
-       LC_SI,
-       0,
-       0,
-       0,
-       0,
-};
-
 /*
  * the moldyn functions
  */
@@ -202,37 +109,27 @@ int set_pressure(t_moldyn *moldyn,double p_ref) {
 
 int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc) {
 
-       moldyn->pt_scalei&=(^(P_SCALE_MASK));
+       moldyn->pt_scale&=(~(P_SCALE_MASK));
        moldyn->pt_scale|=ptype;
        moldyn->p_tc=ptc;
 
-       printf("[moldyn] p/t scaling:\n");
+       printf("[moldyn] p scaling:\n");
 
        printf("  p: %s",ptype?"yes":"no ");
        if(ptype)
                printf(" | type: %02x | factor: %f",ptype,ptc);
        printf("\n");
 
-       printf("  t: %s",ttype?"yes":"no ");
-       if(ttype)
-               printf(" | type: %02x | factor: %f",ttype,ttc);
-       printf("\n");
-
        return 0;
 }
 
 int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc) {
 
-       moldyn->pt_scalei&=(^(T_SCALE_MASK));
+       moldyn->pt_scale&=(~(T_SCALE_MASK));
        moldyn->pt_scale|=ttype;
        moldyn->t_tc=ttc;
 
-       printf("[moldyn] p/t scaling:\n");
-
-       printf("  p: %s",ptype?"yes":"no ");
-       if(ptype)
-               printf(" | type: %02x | factor: %f",ptype,ptc);
-       printf("\n");
+       printf("[moldyn] t scaling:\n");
 
        printf("  t: %s",ttype?"yes":"no ");
        if(ttype)
@@ -323,7 +220,7 @@ int set_potential(t_moldyn *moldyn,u8 type) {
                        moldyn->func3b_k1=tersoff_mult_3bp_k1;
                        moldyn->func3b_j2=tersoff_mult_3bp_j2;
                        moldyn->func3b_k2=tersoff_mult_3bp_k2;
-                       // missing: check 2b bond func
+                       moldyn->check_2b_bond=tersoff_mult_check_2b_bond;
                        break;
                case MOLDYN_POTENTIAL_AM:
                        moldyn->func3b_j1=albe_mult_3bp_j1;
@@ -582,11 +479,16 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
        t_3dvec orig;
        void *ptr;
        t_atom *atom;
+       char name[16];
 
        new=a*b*c;
        count=moldyn->count;
 
        /* how many atoms do we expect */
+       if(type==NONE) {
+               new*=1;
+               printf("[moldyn] WARNING: create 'none' lattice called");
+       }
        if(type==CUBIC) new*=1;
        if(type==FCC) new*=4;
        if(type==DIAMOND) new*=8;
@@ -616,18 +518,21 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
                case CUBIC:
                        set_nn_dist(moldyn,lc);
                        ret=cubic_init(a,b,c,lc,atom,&orig);
+                       strcpy(name,"cubic");
                        break;
                case FCC:
                        if(!origin)
                                v3_scale(&orig,&orig,0.5);
                        set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
                        ret=fcc_init(a,b,c,lc,atom,&orig);
+                       strcpy(name,"fcc");
                        break;
                case DIAMOND:
                        if(!origin)
                                v3_scale(&orig,&orig,0.25);
                        set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
                        ret=diamond_init(a,b,c,lc,atom,&orig);
+                       strcpy(name,"diamond");
                        break;
                default:
                        printf("unknown lattice type (%02x)\n",type);
@@ -644,7 +549,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
        }
 
        moldyn->count+=new;
-       printf("[moldyn] created lattice with %d atoms\n",new);
+       printf("[moldyn] created %s lattice with %d atoms\n",name,new);
 
        for(ret=0;ret<new;ret++) {
                atom[ret].element=element;
@@ -1673,13 +1578,16 @@ int moldyn_integrate(t_moldyn *moldyn) {
                printf("[moldyn] WARNING: cutoff > 0.5 x dim.y\n");
        if(moldyn->cutoff>0.5*moldyn->dim.z)
                printf("[moldyn] WARNING: cutoff > 0.5 x dim.z\n");
-       ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
-       if(ds>0.05*moldyn->nnd)
+       if(moldyn->count) {
+               ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
+               if(ds>0.05*moldyn->nnd)
                printf("[moldyn] WARNING: forces too high / tau too small!\n");
+       }
 
        /* zero absolute time */
-       moldyn->time=0.0;
-       moldyn->total_steps=0;
+       // should have right values!
+       //moldyn->time=0.0;
+       //moldyn->total_steps=0;
 
        /* debugging, ignore */
        moldyn->debug=0;
@@ -1711,7 +1619,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
                temperature_calc(moldyn);
                virial_sum(moldyn);
                pressure_calc(moldyn);
-               //thermodynamic_pressure_calc(moldyn);
+               /*
+               thermodynamic_pressure_calc(moldyn);
+               printf("\n\nDEBUG: numeric pressure calc: %f\n\n",
+                      moldyn->tp/BAR);
+               */
 
                /* calculate fluctuations + averages */
                average_and_fluctuation_calc(moldyn);
@@ -1792,6 +1704,7 @@ printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)",
        sched->count,i,moldyn->total_steps,
        moldyn->t,moldyn->t_avg,
        moldyn->p/BAR,moldyn->p_avg/BAR,
+       //moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
        moldyn->volume,
        (int)(t2.tv_sec-t1.tv_sec));
 
@@ -2422,9 +2335,6 @@ int process_2b_bonds(t_moldyn *moldyn,void *data,
        t_list *this;
 
        lc=&(moldyn->lc);
-
-       link_cell_init(moldyn,VERBOSE);
-
        itom=moldyn->atom;
        
        for(i=0;i<moldyn->count;i++) {
@@ -2684,8 +2594,6 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
        if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
                return 0;
 
-       d=sqrt(d);
-
        /* now count this bonding ... */
        ba=data;