#include "moldyn.h"
#include "report/report.h"
+/* potential includes */
+#include "potentials/harmonic_oscillator.h"
+#include "potentials/lennard_jones.h"
+#include "potentials/albe.h"
+#ifdef TERSOFF_ORIG
+#include "potentials/tersoff_orig.h"
+#else
+#include "potentials/tersoff.h"
+#endif
+
+
/*
* global variables, pse and atom colors (only needed here)
*/
return 0;
}
-int set_potential1b(t_moldyn *moldyn,pf_func1b func) {
-
- moldyn->func1b=func;
-
- return 0;
-}
-
-int set_potential2b(t_moldyn *moldyn,pf_func2b func) {
-
- moldyn->func2b=func;
-
- return 0;
-}
-
-int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func) {
-
- moldyn->func3b_j1=func;
-
- return 0;
-}
-
-int set_potential3b_j2(t_moldyn *moldyn,pf_func2b func) {
-
- moldyn->func3b_j2=func;
-
- return 0;
-}
-
-int set_potential3b_j3(t_moldyn *moldyn,pf_func2b func) {
-
- moldyn->func3b_j3=func;
-
- return 0;
-}
-
-int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func) {
-
- moldyn->func3b_k1=func;
-
- return 0;
-}
-
-int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func) {
-
- moldyn->func3b_k2=func;
-
- return 0;
-}
-
-int set_potential_params(t_moldyn *moldyn,void *params) {
+int set_potential(t_moldyn *moldyn,u8 type) {
- moldyn->pot_params=params;
+ switch(type) {
+ case MOLDYN_POTENTIAL_TM:
+ 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;
+ moldyn->func3b_k2=tersoff_mult_3bp_k2;
+ // missing: check 2b bond func
+ break;
+ case MOLDYN_POTENTIAL_AM:
+ moldyn->func3b_j1=albe_mult_3bp_j1;
+ moldyn->func3b_k1=albe_mult_3bp_k1;
+ moldyn->func3b_j2=albe_mult_3bp_j2;
+ moldyn->func3b_k2=albe_mult_3bp_k2;
+ moldyn->check_2b_bond=albe_mult_check_2b_bond;
+ break;
+ case MOLDYN_POTENTIAL_HO:
+ moldyn->func2b=harmonic_oscillator;
+ moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond;
+ break;
+ case MOLDYN_POTENTIAL_LJ:
+ moldyn->func2b=lennard_jones;
+ moldyn->check_2b_bond=lennard_jones_check_2b_bond;
+ break;
+ default:
+ printf("[moldyn] set potential: unknown type %02x\n",
+ type);
+ return -1;
+ }
return 0;
}
return 0;
}
+int moldyn_free_save_file(t_moldyn *moldyn) {
+
+ free(moldyn->atom);
+
+ return 0;
+}
+
int moldyn_load(t_moldyn *moldyn) {
// later ...
return 0;
}
+/*
+ * function to find/callback all combinations of 2 body bonds
+ */
+
+int process_2b_bonds(t_moldyn *moldyn,void *data,
+ int (*process)(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
+ void *data,u8 bc)) {
+
+ t_linkcell *lc;
+#ifdef STATIC_LISTS
+ int *neighbour[27];
+ int p;
+#else
+ t_list neighbour[27];
+#endif
+ u8 bc;
+ t_atom *itom,*jtom;
+ int i,j;
+ t_list *this;
+
+ lc=&(moldyn->lc);
+
+ link_cell_init(moldyn,VERBOSE);
+
+ itom=moldyn->atom;
+
+ for(i=0;i<moldyn->count;i++) {
+ /* neighbour indexing */
+ link_cell_neighbour_index(moldyn,
+ (itom[i].r.x+moldyn->dim.x/2)/lc->x,
+ (itom[i].r.y+moldyn->dim.y/2)/lc->x,
+ (itom[i].r.z+moldyn->dim.z/2)/lc->x,
+ neighbour);
+
+ for(j=0;j<27;j++) {
+
+ bc=(j<lc->dnlc)?0:1;
+
+#ifdef STATIC_LISTS
+ p=0;
+
+ while(neighbour[j][p]!=0) {
+
+ jtom=&(moldyn->atom[neighbour[j][p]]);
+ p++;
+#else
+ this=&(neighbour[j]);
+ list_reset_f(this);
+
+ if(this->start==NULL)
+ continue;
+
+ do {
+
+ jtom=this->current->data;
+#endif
+
+ /* process bond */
+ process(moldyn,&(itom[i]),jtom,data,bc);
+
+#ifdef STATIC_LISTS
+ }
+#else
+ } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+#endif
+ }
+ }
+
+ return 0;
+
+}
+
/*
* post processing functions
*/
return 0;
}
-int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
+int bonding_analyze(t_moldyn *moldyn,double *cnt) {
+
+ return 0;
+}
+
+int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom,
+ t_atom *jtom,void *data,u8 bc) {
- int slots;
- double *stat;
- int i,j;
- t_linkcell *lc;
-#ifdef STATIC_LISTS
- int *neighbour[27];
- int p;
-#else
- t_list neighbour[27];
-#endif
- t_atom *itom,*jtom;
- t_list *this;
- unsigned char bc;
t_3dvec dist;
double d;
- double norm;
- int o,s;
- unsigned char ibrand;
+ int s;
+ t_pcc *pcc;
- lc=&(moldyn->lc);
+ /* only count pairs once,
+ * skip same atoms */
+ if(itom->tag>=jtom->tag)
+ return 0;
+
+ /*
+ * pair correlation calc
+ */
+
+ /* get pcc data */
+ pcc=data;
+
+ /* distance */
+ v3_sub(&dist,&(jtom->r),&(itom->r));
+ if(bc) check_per_bound(moldyn,&dist);
+ d=v3_absolute_square(&dist);
+
+ /* ignore if greater cutoff */
+ if(d>moldyn->cutoff_square)
+ return 0;
+
+ /* fill the slots */
+ d=sqrt(d);
+ s=(int)(d/pcc->dr);
+
+ /* should never happen but it does 8) -
+ * related to -ffloat-store problem! */
+ if(s>=pcc->o1) {
+ printf("[moldyn] WARNING: pcc (%d/%d)",
+ s,pcc->o1);
+ printf("\n");
+ s=pcc->o1-1;
+ }
+
+ if(itom->brand!=jtom->brand) {
+ /* mixed */
+ pcc->stat[s]+=1;
+ }
+ else {
+ /* type a - type a bonds */
+ if(itom->brand==0)
+ pcc->stat[s+pcc->o1]+=1;
+ else
+ /* type b - type b bonds */
+ pcc->stat[s+pcc->o2]+=1;
+ }
+
+ return 0;
+}
+
+int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
+
+ t_pcc pcc;
+ double norm;
+ int i;
- slots=2.0*moldyn->cutoff/dr;
- o=2*slots;
+ pcc.dr=dr;
+ pcc.o1=moldyn->cutoff/dr;
+ pcc.o2=2*pcc.o1;
- if(slots*dr<=moldyn->cutoff)
+ if(pcc.o1*dr<=moldyn->cutoff)
printf("[moldyn] WARNING: pcc (low #slots)\n");
printf("[moldyn] pair correlation calc info:\n");
printf(" temperature: cur=%f avg=%f\n",moldyn->t,moldyn->t_avg);
if(ptr!=NULL) {
- stat=(double *)ptr;
+ pcc.stat=(double *)ptr;
}
else {
- stat=(double *)malloc(3*slots*sizeof(double));
- if(stat==NULL) {
+ pcc.stat=(double *)malloc(3*pcc.o1*sizeof(double));
+ if(pcc.stat==NULL) {
perror("[moldyn] pair correlation malloc");
return -1;
}
}
- memset(stat,0,3*slots*sizeof(double));
+ memset(pcc.stat,0,3*pcc.o1*sizeof(double));
- link_cell_init(moldyn,VERBOSE);
+ /* process */
+ process_2b_bonds(moldyn,&pcc,calculate_pair_correlation_process);
- itom=moldyn->atom;
-
- for(i=0;i<moldyn->count;i++) {
- /* neighbour indexing */
- link_cell_neighbour_index(moldyn,
- (itom[i].r.x+moldyn->dim.x/2)/lc->x,
- (itom[i].r.y+moldyn->dim.y/2)/lc->x,
- (itom[i].r.z+moldyn->dim.z/2)/lc->x,
- neighbour);
+ /* normalization */
+ for(i=1;i<pcc.o1;i++) {
+ // normalization: 4 pi r^2 dr
+ // here: not double counting pairs -> 2 pi r r dr
+ // ... and actually it's a constant times r^2
+ norm=i*i*dr*dr;
+ pcc.stat[i]/=norm;
+ pcc.stat[pcc.o1+i]/=norm;
+ pcc.stat[pcc.o2+i]/=norm;
+ }
+ /* */
- /* brand of atom i */
- ibrand=itom[i].brand;
-
- for(j=0;j<27;j++) {
+ if(ptr==NULL) {
+ /* todo: store/print pair correlation function */
+ free(pcc.stat);
+ }
- bc=(j<lc->dnlc)?0:1;
+ return 0;
+}
-#ifdef STATIC_LISTS
- p=0;
+int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
+ void *data,u8 bc) {
- while(neighbour[j][p]!=0) {
+ t_ba *ba;
+ t_3dvec dist;
+ double d;
- jtom=&(moldyn->atom[neighbour[j][p]]);
- p++;
-#else
- this=&(neighbour[j]);
- list_reset_f(this);
+ if(itom->tag>=jtom->tag)
+ return 0;
- if(this->start==NULL)
- continue;
+ /* distance */
+ v3_sub(&dist,&(jtom->r),&(itom->r));
+ if(bc) check_per_bound(moldyn,&dist);
+ d=v3_absolute_square(&dist);
- do {
+ /* ignore if greater or equal cutoff */
+ if(d>moldyn->cutoff_square)
+ return 0;
- jtom=this->current->data;
-#endif
- /* only count pairs once,
- * skip same atoms */
- if(itom[i].tag>=jtom->tag)
- continue;
+ /* check for potential bond */
+ if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
+ return 0;
- /*
- * pair correlation calc
- */
+ d=sqrt(d);
- /* distance */
- v3_sub(&dist,&(jtom->r),&(itom[i].r));
- if(bc) check_per_bound(moldyn,&dist);
- d=v3_absolute_square(&dist);
+ /* now count this bonding ... */
+ ba=data;
- /* ignore if greater or equal cutoff */
- if(d>=4.0*moldyn->cutoff_square)
- continue;
+ /* increase total bond counter
+ * ... double counting!
+ */
+ ba->tcnt+=2;
- /* fill the slots */
- d=sqrt(d);
- s=(int)(d/dr);
-
- /* should never happen but it does 8) -
- * related to -ffloat-store problem! */
- if(s>=slots) {
- printf("[moldyn] WARNING: pcc (%d/%d)",
- s,slots);
- printf("\n");
- s=slots-1;
- }
+ if(itom->brand==0)
+ ba->acnt[jtom->tag]+=1;
+ else
+ ba->bcnt[jtom->tag]+=1;
+
+ if(jtom->brand==0)
+ ba->acnt[itom->tag]+=1;
+ else
+ ba->bcnt[itom->tag]+=1;
- if(ibrand!=jtom->brand) {
- /* mixed */
- stat[s]+=1;
- }
- else {
- /* type a - type a bonds */
- if(ibrand==0)
- stat[s+slots]+=1;
- else
- /* type b - type b bonds */
- stat[s+o]+=1;
- }
-#ifdef STATIC_LISTS
- }
-#else
- } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
-#endif
- }
- }
+ return 0;
+}
- /* normalization */
- for(i=1;i<slots;i++) {
- // normalization: 4 pi r^2 dr
- // here: not double counting pairs -> 2 pi r r dr
- // ... and actually it's a constant times r^2
- norm=i*i*dr*dr;
- stat[i]/=norm;
- stat[slots+i]/=norm;
- stat[o+i]/=norm;
+int bond_analyze(t_moldyn *moldyn,double *quality) {
+
+ // by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total
+
+ int qcnt;
+ int ccnt,cset;
+ t_ba ba;
+ int i;
+ t_atom *atom;
+
+ ba.acnt=malloc(moldyn->count*sizeof(int));
+ if(ba.acnt==NULL) {
+ perror("[moldyn] bond analyze malloc (a)");
+ return -1;
}
- /* */
+ memset(ba.acnt,0,moldyn->count*sizeof(int));
- if(ptr==NULL) {
- /* todo: store/print pair correlation function */
- free(stat);
+ ba.bcnt=malloc(moldyn->count*sizeof(int));
+ if(ba.bcnt==NULL) {
+ perror("[moldyn] bond analyze malloc (b)");
+ return -1;
}
+ memset(ba.bcnt,0,moldyn->count*sizeof(int));
- free(moldyn->atom);
+ ba.tcnt=0;
+ qcnt=0;
+ ccnt=0;
+ cset=0;
- link_cell_shutdown(moldyn);
+ atom=moldyn->atom;
- return 0;
-}
+ process_2b_bonds(moldyn,&ba,bond_analyze_process);
-int analyze_bonds(t_moldyn *moldyn) {
+ for(i=0;i<moldyn->count;i++) {
+ if(atom[i].brand==0) {
+ if((ba.acnt[i]==0)&(ba.bcnt[i]==4))
+ qcnt+=4;
+ }
+ else {
+ if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) {
+ qcnt+=4;
+ ccnt+=1;
+ }
+ cset+=1;
+ }
+ }
-
-
+ printf("[moldyn] bond analyze: c_cnt=%d | set=%d\n",ccnt,cset);
+ printf("[moldyn] bond analyze: q_cnt=%d | tot=%d\n",qcnt,ba.tcnt);
+
+ if(quality) {
+ quality[0]=1.0*ccnt/cset;
+ quality[1]=1.0*qcnt/ba.tcnt;
+ }
+ else {
+ printf("[moldyn] bond analyze: c_bnd_q=%f\n",1.0*qcnt/ba.tcnt);
+ printf("[moldyn] bond analyze: tot_q=%f\n",1.0*qcnt/ba.tcnt);
+ }
return 0;
}