#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)
*/
int set_cutoff(t_moldyn *moldyn,double cutoff) {
moldyn->cutoff=cutoff;
+ moldyn->cutoff_square=cutoff*cutoff;
printf("[moldyn] cutoff [A]: %f\n",moldyn->cutoff);
return 0;
}
-int set_bondlen(t_moldyn *moldyn,double b0,double b1,double bm) {
-
- moldyn->bondlen[0]=b0*b0;
- moldyn->bondlen[1]=b1*b1;
- if(bm<0)
- moldyn->bondlen[2]=b0*b1;
- else
- moldyn->bondlen[2]=bm*bm;
-
- return 0;
-}
-
int set_temperature(t_moldyn *moldyn,double t_ref) {
moldyn->t_ref=t_ref;
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;
+int set_potential(t_moldyn *moldyn,u8 type) {
- 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) {
-
- 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;
}
/* sqaure of some variables */
moldyn->tau_square=moldyn->tau*moldyn->tau;
- moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
/* get current time */
gettimeofday(&t1,NULL);
ret=read(fd,line+count,1);
if(ret<=0) return ret;
if(line[count]=='\n') {
- line[count]='\0';
+ memset(line+count,0,max-count-1);
+ //line[count]='\0';
return count+1;
}
count+=1;
if(bc) check_per_bound(moldyn,&dist);
d=v3_absolute_square(&dist);
- /* ignore if greater or equal 2 times cutoff */
- if(d>=4.0*moldyn->cutoff_square)
+ /* ignore if greater cutoff */
+ if(d>moldyn->cutoff_square)
return 0;
/* fill the slots */
int i;
pcc.dr=dr;
- pcc.o1=2.0*moldyn->cutoff/dr;
+ pcc.o1=moldyn->cutoff/dr;
pcc.o2=2*pcc.o1;
if(pcc.o1*dr<=moldyn->cutoff)
d=v3_absolute_square(&dist);
/* ignore if greater or equal cutoff */
- if(d>=moldyn->cutoff_square)
+ if(d>moldyn->cutoff_square)
return 0;
+ /* check for potential bond */
+ if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
+ return 0;
+
+ d=sqrt(d);
+
/* now count this bonding ... */
ba=data;
// 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;
memset(ba.bcnt,0,moldyn->count*sizeof(int));
ba.tcnt=0;
-
qcnt=0;
+ ccnt=0;
+ cset=0;
atom=moldyn->atom;
qcnt+=4;
}
else {
- if((ba.acnt[i]==4)&(ba.bcnt[i]==0))
+ if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) {
qcnt+=4;
+ ccnt+=1;
+ }
+ cset+=1;
}
}
-printf("%d %d\n",qcnt,ba.tcnt);
- if(quality)
- *quality=1.0*qcnt/ba.tcnt;
- else
- printf("[moldyn] bond analyze: quality = %f\n",
- 1.0*qcnt/ba.tcnt);
+ 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;
}
return 0;
}
+int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
+ void *data,u8 bc) {
+
+ t_vb *vb;
+
+ vb=data;
+
+ if(itom->tag>=jtom->tag)
+ return 0;
+
+ if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
+ return 0;
+
+ if((itom->attr&ATOM_ATTR_VB)|(jtom->attr&ATOM_ATTR_VB))
+ dprintf(vb->fd,"# [B] %f %f %f %f %f %f\n",
+ itom->r.x,itom->r.y,itom->r.z,
+ jtom->r.x,jtom->r.y,jtom->r.z);
+
+ return 0;
+}
+
int visual_atoms(t_moldyn *moldyn) {
- int i,j,fd;
+ int i;
char file[128+64];
t_3dvec dim;
double help;
t_visual *v;
t_atom *atom;
- t_atom *btom;
- t_linkcell *lc;
-#ifdef STATIC_LISTS
- int *neighbour[27];
- int p;
-#else
- t_list neighbour[27];
-#endif
- u8 bc;
- t_3dvec dist;
- double d2;
- u8 brand;
+ t_vb vb;
v=&(moldyn->vis);
dim.x=v->dim.x;
dim.y=v->dim.y;
dim.z=v->dim.z;
atom=moldyn->atom;
- lc=&(moldyn->lc);
help=(dim.x+dim.y);
sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time);
- fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
- if(fd<0) {
+ vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+ if(vb.fd<0) {
perror("open visual save file fd");
return -1;
}
/* write the actual data file */
// povray header
- dprintf(fd,"# [P] %d %07.f <%f,%f,%f>\n",
+ dprintf(vb.fd,"# [P] %d %07.f <%f,%f,%f>\n",
moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help);
// atomic configuration
- for(i=0;i<moldyn->count;i++) {
+ for(i=0;i<moldyn->count;i++)
// atom type, positions, color and kinetic energy
- dprintf(fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element],
- atom[i].r.x,
- atom[i].r.y,
- atom[i].r.z,
- pse_col[atom[i].element],
- atom[i].ekin);
-
- /*
- * bond detection should usually be done by potential
- * functions. brrrrr! EVIL!
- *
- * todo: potentials need to export a 'find_bonds' function!
- */
-
- // bonds between atoms
- if(!(atom[i].attr&ATOM_ATTR_VB))
- continue;
- link_cell_neighbour_index(moldyn,
- (atom[i].r.x+moldyn->dim.x/2)/lc->x,
- (atom[i].r.y+moldyn->dim.y/2)/lc->y,
- (atom[i].r.z+moldyn->dim.z/2)/lc->z,
- neighbour);
- for(j=0;j<27;j++) {
- bc=j<lc->dnlc?0:1;
-#ifdef STATIC_LISTS
- p=0;
- while(neighbour[j][p]!=0) {
- btom=&(atom[neighbour[j][p]]);
- p++;
-#else
- list_reset_f(&neighbour[j]);
- if(neighbour[j].start==NULL)
- continue;
- do {
- btom=neighbour[j].current->data;
-#endif
- if(btom==&atom[i]) // skip identical atoms
- continue;
- //if(btom<&atom[i]) // skip half of them
- // continue;
- v3_sub(&dist,&(atom[i].r),&(btom->r));
- if(bc) check_per_bound(moldyn,&dist);
- d2=v3_absolute_square(&dist);
- brand=atom[i].brand;
- if(brand==btom->brand) {
- if(d2>moldyn->bondlen[brand])
- continue;
- }
- else {
- if(d2>moldyn->bondlen[2])
- continue;
- }
- dprintf(fd,"# [B] %f %f %f %f %f %f\n",
- atom[i].r.x,atom[i].r.y,atom[i].r.z,
- btom->r.x,btom->r.y,btom->r.z);
-#ifdef STATIC_LISTS
- }
-#else
- } while(list_next_f(&neighbour[j])!=L_NO_NEXT_ELEMENT);
-#endif
- }
- }
-
+ dprintf(vb.fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element],
+ atom[i].r.x,
+ atom[i].r.y,
+ atom[i].r.z,
+ pse_col[atom[i].element],
+ atom[i].ekin);
+
+ // bonds between atoms
+ process_2b_bonds(moldyn,&vb,visual_bonds_process);
+
// boundaries
if(dim.x) {
- dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
-dim.x/2,-dim.y/2,-dim.z/2,
dim.x/2,-dim.y/2,-dim.z/2);
- dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
-dim.x/2,-dim.y/2,-dim.z/2,
-dim.x/2,dim.y/2,-dim.z/2);
- dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
dim.x/2,dim.y/2,-dim.z/2,
dim.x/2,-dim.y/2,-dim.z/2);
- dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
-dim.x/2,dim.y/2,-dim.z/2,
dim.x/2,dim.y/2,-dim.z/2);
- dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
-dim.x/2,-dim.y/2,dim.z/2,
dim.x/2,-dim.y/2,dim.z/2);
- dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
-dim.x/2,-dim.y/2,dim.z/2,
-dim.x/2,dim.y/2,dim.z/2);
- dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
dim.x/2,dim.y/2,dim.z/2,
dim.x/2,-dim.y/2,dim.z/2);
- dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
-dim.x/2,dim.y/2,dim.z/2,
dim.x/2,dim.y/2,dim.z/2);
- dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
-dim.x/2,-dim.y/2,dim.z/2,
-dim.x/2,-dim.y/2,-dim.z/2);
- dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
-dim.x/2,dim.y/2,dim.z/2,
-dim.x/2,dim.y/2,-dim.z/2);
- dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
dim.x/2,-dim.y/2,dim.z/2,
dim.x/2,-dim.y/2,-dim.z/2);
- dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
dim.x/2,dim.y/2,dim.z/2,
dim.x/2,dim.y/2,-dim.z/2);
}
- close(fd);
+ close(vb.fd);
return 0;
}
projects / physik / posic.git / blobdiff