int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
int i;
- t_ho_params hop;
- t_lj_params ljp;
- t_tersoff_params tp;
- double s,e;
memset(moldyn,0,sizeof(t_moldyn));
moldyn->tau=MOLDYN_TAU;
moldyn->time_steps=MOLDYN_RUNS;
moldyn->integrate=velocity_verlet;
- moldyn->potential_force_function=lennard_jones;
/* parse argv */
for(i=1;i<argc;i++) {
int moldyn_integrate(t_moldyn *moldyn) {
- int i;
+ int i,sched;
unsigned int e,m,s,d,v;
t_3dvec p;
/* calculate initial forces */
moldyn->potential_force_function(moldyn);
+ for(sched=0;sched<moldyn->schedule.content_count;sched++) {
+ moldyn->tau=;
+ moldyn->tau_square=;
+
+ // hier weiter ...
+
for(i=0;i<moldyn->time_steps;i++) {
/* integration step */
* tersoff potential & force for 2 sorts of atoms
*/
-/* tersoff 2 body part */
+/* tersoff 1 body part */
+int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
+
+ int num;
+ t_tersoff_mult_params *params;
+ t_tersoff_exchange *exchange;
+
+ num=ai->bnum;
+ params=moldyn->pot1b_params;
+ exchange=&(params->exchange);
+
+ /*
+ * simple: point constant parameters only depending on atom i to
+ * their right values
+ */
-int tersoff_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+ exchange->beta=&(params->beta[num]);
+ exchange->n=&(params->n[num]);
+ exchange->c=&(params->c[num]);
+ exchange->d=&(params->d[num]);
+ exchange->h=&(params->h[num]);
- t_tersoff_params *params;
+ exchange->betan=pow(*(exchange->beta),*(exchange->n));
+ exchange->c2=params->c[num]*params->c[num];
+ exchange->d2=params->d[num]*params->d[num];
+ exchange->c2d2=exchange->c2/exchange->d2;
+
+ return 0;
+}
+
+/* tersoff 2 body part */
+int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+
+ t_tersoff_mult_params *params;
+ t_tersoff_exchange *exchange;
t_3dvec dist_ij;
double d_ij;
+ double A,B,R,S,lambda;
+ int num;
params=moldyn->pot_params;
+ num=ai->bnum;
+ exchange=&(params->exchange);
+
+ exchange->run3bp=0;
/*
* we need: f_c, df_c, f_r, df_r
*
- * therefore we need: R, S
+ * therefore we need: R, S, A, lambda
*/
v3_sub(&dist_ij,&(ai->r),&(aj->r));
if(bc) check_per_bound(moldyn,&dist_ij);
- if(ai->bnum==aj->bnum) {
- S=params->S[ai->bnum];
- R=params->R[ai->bnum];
+ /* save for use in 3bp */ /* REALLY ?!?!?! */
+ exchange->dist_ij=dist_ij;
+
+ /* constants */
+ if(num==aj->bnum) {
+ S=params->S[num];
+ R=params->R[num];
+ A=params->A[num];
+ lambda=params->lambda[num];
+ /* more constants depending of atoms i and j, needed in 3bp */
+ params->exchange.B=&(params->B[num]);
+ params->exchange.mu=params->mu[num];
+ params->exchange.chi=1.0;
}
else {
S=params->Smixed;
R=params->Rmixed;
+ A=params->Amixed;
+ lambda=params->lambda_m;
+ /* more constants depending of atoms i and j, needed in 3bp */
+ params->exchange.B=&(params->Bmixed);
+ params->exchange.mu=&(params->mu_m);
+ params->exchange.chi=params->chi;
}
d_ij=v3_norm(&dist_ij);
- if(d_ij<=S) {
- f_r=A*exp(-lamda*d_ij);
- df_r=-lambda*f_r/d_ij;
- if(d_ij<R) {
- /* f_c = 1, df_c = 0 */
- v3_scale(&force,&dist_ij,df_r);
- }
- else {
- s_r=S-R;
- arg=PI*(d_ij-R)/s_r;
- f_c=0.5+0.5*cos(arg);
- df_c=-0.5*sin(arg)*(PI/(s_r*d_ij));
- scale=df_c*f_r+df_r*f_c;
- v3_scale(&force,&dist_ij,scale);
- v3_add(&(ai->f),&(ai->f),&force);
- }
- moldyn->energy+=(f_r*f_c);
- }
-
- return 0;
-}
-
-/* tersoff 3 body part */
-
-int tersoff(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc,u8 bck) {
-
- t_tersoff_params *params;
- t_3dvec dist_ij;
- double d_ij;
-
- params=moldyn->pot_params;
-
- /* 2 body part of the tersoff potential */
-
- v3_sub(&dist_ij,&(ai->r),&(aj->r));
- if(bc) check_per_bound(moldyn,&dist_ij);
- d_ij=v3_norm(&dist_ij);
- if(d_ij<=S) {
-
- /* determine the tersoff parameters */
- if(atom[i].element!=btom->element) {
- S=sqrt(TERSOFF_S[e1]*TERSOFF_S[e2]);
- R=R_m;
- A=;
- lambda=;
- B=;
- mu=;
- chi=;
- beta=;
- betaN=;
-
- if(d_ij<=R) {
- df_r=-lambda*A*exp(-lambda*d_ij)/d_ij;
- v3_scale(&force,&dist_ij,df_r);
- v3_add(&(atom[i].f),&(atom[i].f),
- &force);
- }
- else {
- s_r=S-R;
- arg1=PI*(d_ij-R)/s_r;
- f_c=0.5+0.5*cos(arg1);
- df_c=-0.5*sin(arg1)*(PI/(s_r*d_ij));
- f_r=A*exp(-lambda*d_ij);
- df_r=-lambda*f_r/d_ij;
- scale=df_c*f_r+df_r*f_c;
- v3_scale(&force,&dist_ij,scale);
- v3_add(&(atom[i].f),&(atom[i].f),
- &force);
- }
- }
- else
- continue;
-
-
- /* end 2 body stuff */
-
- /* determine cell neighbours of btom */
- ki=(btom->r.x+(moldyn->dim.x/2))/lc->x;
- kj=(btom->r.y+(moldyn->dim.y/2))/lc->y;
- kk=(btom->r.z+(moldyn->dim.z/2))/lc->z;
- ck=link_cell_neighbour_index(moldyn,ki,kj,kk,
- neighbourk);
-
- /* go for zeta - 3 body stuff! */
- zeta=0.0;
- d_ij2=d_ij*d_ij;
-
- /* cell of btom */
- thisk=&(neighbourk[0]);
- list_reset(thisk);
- do {
- ktom=thisk->current->data;
- if(ktom==btom)
- continue;
- if(ktom==&(atom[i]))
- continue;
-
- /* 3 body stuff (1) */
-
- v3_sub(&dist_ik,ktom,&(atom[i]));
- d_ik=v3_norm(&dist_ik);
- if(d_ik<=Sik) {
-
- Rik=;
- Sik=;
- Aik=;
- lambda_ik=;
- Bik=;
- mu_ik=;
- omega_ik=;
- c_i=;
- d_i=;
- h_i=;
-
-
- if(d_ik<=Rik) {
- f_cik=1.0;
- df_cik=0.0;
- }
- else {
- sik_rik=Sik-Rik;
- arg1ik=PI*(d_ik-Rik)/sik_rik;
- f_cik=0.5+0.5*cos(arg1ik);
- df_cik=-0.5*sin(arg1ik)* \
- (PI/(sik_rik*d_ik));
- f_rik=Aik*exp(-lambda_ik*d_ik);
- f_aik=-Bik*exp(-mu_ik*d_ik);
- }
-
- v3_sub(&distance_jk,ktom,btom);
- cos_theta=(d_ij2+d_ik*d_ik-d_jk*d_jk)/\
- (2*d_ij*d_ik);
- sin_theta=sqrt(1.0/\
- (cos_theta*cos_theta));
- theta=arccos(cos_theta);
-
-
- }
- else
- continue;
-
- /* end 3 body stuff (1) */
-
-
- } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
-
- /* direct neighbours of btom cell */
- for(k=1;k<ck;k++) {
- thisk=&(neighbourk[k]);
- list_reset(thisk);
- if(thisk->start!=NULL) {
-
- do {
- ktom=thisk->current->data;
- if(ktom==&(atom[i]))
- continue;
-
- /* 3 body stuff (2) */
-
- } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
-
- }
- }
-
- /* indirect neighbours of btom cell */
- for(k=ck;k<27;k++) {
- thisk=&(neighbourk[k]);
- list_reset(thisk);
- if(thisk->start!=NULL) {
-
- do {
- ktom=thisk->current->data;
- if(ktom==&(atom[i]))
- continue;
-
- /* 3 body stuff */
-
- } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
-
- }
- }
-
-
- } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+ /* save for use in 3bp */
+ exchange->d_ij=d_ij;
- /*
- * direct neighbour cells of atom i
- */
- for(j=1;j<c;j++) {
- this=&(neighbour[j]);
- list_reset(this);
- if(this->start!=NULL) {
+ if(d_ij>S)
+ return 0;
- do {
- btom=this->current->data;
+ f_r=A*exp(-lamda*d_ij);
+ df_r=-lambda*f_r/d_ij;
- /* 2 body stuff */
+ /* f_a, df_a calc + save for 3bp use */
+ exchange->f_a=-B*exp(-mu*d_ij);
+ exchange->df_a=-mu*exchange->f_a/d_ij;
+ if(d_ij<R) {
+ /* f_c = 1, df_c = 0 */
+ f_c=1.0;
+ df_c=0.0;
+ v3_scale(&force,&dist_ij,df_r);
+ }
+ else {
+ s_r=S-R;
+ arg=PI*(d_ij-R)/s_r;
+ f_c=0.5+0.5*cos(arg);
+ df_c=-0.5*sin(arg)*(PI/(s_r*d_ij));
+ scale=df_c*f_r+df_r*f_c;
+ v3_scale(&force,&dist_ij,scale);
+ }
- /* determine cell neighbours of btom */
- ki=(btom->r.x+(moldyn->dim.x/2))/lc->x;
- kj=(btom->r.y+(moldyn->dim.y/2))/lc->y;
- kk=(btom->r.z+(moldyn->dim.z/2))/lc->z;
- ck=link_cell_neighbour_index(moldyn,ki,kj,kk,
- neighbourk);
-
- /* cell of btom */
- thisk=&(neighbourk[0]);
- list_reset(thisk);
- do {
- ktom=thisk->current->data;
- if(ktom==btom)
- continue;
- if(ktom==&(atom[i]))
- continue;
-
- /* 3 body stuff (1) */
-
- } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
-
- /* direct neighbours of btom cell */
- for(k=1;k<ck;k++) {
- thisk=&(neighbourk[k]);
- list_reset(thisk);
- if(thisk->start!=NULL) {
-
- do {
- ktom=thisk->current->data;
- if(ktom==&(atom[i]))
- continue;
-
- /* 3 body stuff (2) */
-
- } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
-
- }
- }
-
- /* indirect neighbours of btom cell */
- for(k=ck;k<27;k++) {
- thisk=&(neighbourk[k]);
- list_reset(thisk);
- if(thisk->start!=NULL) {
-
- do {
- ktom=thisk->current->data;
- if(ktom==&(atom[i]))
- continue;
-
- /* 3 body stuff (3) */
-
- } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
-
- }
- }
-
-
- } while(list_next(this)!=L_NO_NEXT_ELEMENT);
-
- }
- }
-
- /*
- * indirect neighbour cells of atom i
- */
- for(j=c;j<27;j++) {
- this=&(neighbour[j]);
- list_reset(this);
- if(this->start!=NULL) {
-
- do {
- btom=this->current->data;
-
- /* 2 body stuff */
-
-
- /* determine cell neighbours of btom */
- ki=(btom->r.x+(moldyn->dim.x/2))/lc->x;
- kj=(btom->r.y+(moldyn->dim.y/2))/lc->y;
- kk=(btom->r.z+(moldyn->dim.z/2))/lc->z;
- ck=link_cell_neighbour_index(moldyn,ki,kj,kk,
- neighbourk);
-
- /* cell of btom */
- thisk=&(neighbourk[0]);
- list_reset(thisk);
- do {
- ktom=thisk->current->data;
- if(ktom==btom)
- continue;
- if(ktom==&(atom[i]))
- continue;
-
- /* 3 body stuff (1) */
-
- } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
+ /* add forces */
+ v3_add(&(ai->f),&(ai->f),&force);
+ /* energy is 0.5 f_r f_c, but we will sum it up twice ... */
+ moldyn->energy+=(0.25*f_r*f_c);
- /* direct neighbours of btom cell */
- for(k=1;k<ck;k++) {
- thisk=&(neighbourk[k]);
- list_reset(thisk);
- if(thisk->start!=NULL) {
+ /* save for use in 3bp */
+ exchange->f_c=f_c;
+ exchange->df_c=df_c;
- do {
- ktom=thisk->current->data;
- if(ktom==&(atom[i]))
- continue;
+ /* enable the run of 3bp function */
+ exchange->run3bp=1;
- /* 3 body stuff (2) */
+ return 0;
+}
- } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
+/* tersoff 3 body part */
- }
- }
+int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
- /* indirect neighbours of btom cell */
- for(k=ck;k<27;k++) {
- thisk=&(neighbourk[k]);
- list_reset(thisk);
- if(thisk->start!=NULL) {
+ t_tersoff_mult_params *params;
+ t_tersoff_exchange *exchange;
+ t_3dvec dist_ij,dist_ik,dist_jk;
+ t_3dvec temp,force;
+ double R,S,s_r;
+ double d_ij,d_ik,d_jk;
+ double f_c,df_c,b_ij,f_a,df_a;
+ double n,c,d,h,neta,betan,betan_1;
+ double theta,cos_theta,sin_theta;
+ int num;
- do {
- ktom=thisk->current->data;
- if(ktom==&(atom[i]))
- continue;
+ params=moldyn->pot_params;
+ num=ai->bnum;
+ exchange=params->exchange;
- /* 3 body stuff (3) */
+ if(!(exchange->run3bp))
+ return 0;
- } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
+ /*
+ * we need: f_c, d_fc, b_ij, db_ij, f_a, df_a
+ *
+ * we got f_c, df_c, f_a, df_a from 2bp calculation
+ */
- }
- }
+ d_ij=exchange->d_ij;
+ d_ij2=exchange->d_ij2;
+ f_a=params->exchange.f_a;
+ df_a=params->exchange.df_a;
+
+ /* d_ij is <= S, as we didn't return so far! */
- } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+ /*
+ * calc of b_ij (scalar) and db_ij (vector)
+ *
+ * - for b_ij: chi, beta, f_c_ik, w(=1), c, d, h, n, cos_theta
+ *
+ * - for db_ij: d_theta, sin_theta, cos_theta, f_c_ik, df_c_ik,
+ * w_ik,
+ *
+ */
- }
- }
-
+
+ v3_sub(&dist_ik,&(aj->i),&(ak->r));
+ if(bc) check_per_bound(moldyn,&dist_ik);
+ d_ik=v3_norm(&dist_ik);
+
+ /* constants for f_c_ik calc */
+ if(num==ak->bnum) {
+ R=params->R[num];
+ S=params->S[num];
+ }
+ else {
+ R=params->Rmixed;
+ S=params->Smixed;
}
- moldyn->energy=0.5*u;
+ /* calc of f_c_ik */
+ if(d_ik>S)
+ return 0;
+ if(d_ik<R) {
+ /* f_c_ik = 1, df_c_ik = 0 */
+ f_c_ik=1.0;
+ df_c_ik=0.0;
+ }
+ else {
+ s_r=S-R;
+ arg=PI*(d_ik-R)/s_r;
+ f_c_ik=0.5+0.5*cos(arg);
+ df_c_ik=-0.5*sin(arg)*(PI/(s_r*d_ik));
+ }
+
+ v3_sub(&dist_jk,&(aj->r),&(ak->r));
+ if(bc) check_per_bound(moldyn,&dist_jk);
+ d_jk=v3_norm(&dist_jk);
+
+ beta=*(exchange->beta);
+ betan=exchange->betan;
+ n=*(exchange->n);
+ c=*(exchange->c);
+ d=*(exchange->d);
+ h=*(exchange->h);
+ c2=exchange->c2;
+ d2=exchange->d2;
+ c2d2=exchange->c2d2;
+
+ numer=d_ij2+d_ik*d_ik-d_jk*d_jk;
+ denom=2*d_ij*d_ik;
+ cos_theta=numer/denom;
+ sin_theta=sqrt(1.0-(cos_theta*cos_theta));
+ theta=arccos(cos_theta);
+ d_theta=(-1.0/sqrt(1.0-cos_theta*cos_theta))/(denom*denom);
+ d_theta1=2*denom-numer*2*d_ik/d_ij;
+ d_theta2=2*denom-numer*2*d_ij/d_ik;
+ d_theta1*=d_theta;
+ d_theta2*=d_theta;
+
+ h_cos=(h-cos_theta);
+ h_cos2=h_cos*h_cos;
+ d2_h_cos2=d2-h_cos2;
+
+ /* some usefull expressions */
+ frac1=c2/(d2-h_cos2);
+ bracket1=1+c2d2-frac1;
+ bracket2=f_c_ik*bracket1;
+ bracket2_n_1=pow(bracket2,n-1.0);
+ bracket2_n=bracket2_n_1*bracket2;
+ bracket3=1+betan*bracket2_n;
+ bracket3_pow_1=pow(bracket3,(-1.0/(2.0*n))-1.0);
+ bracket3_pow=bracket3_pow_1*bracket3;
+
+ /* now go on with calc of b_ij and derivation of b_ij */
+ b_ij=chi*bracket3_pow;
+
+ /* derivation of theta */
+ v3_scale(&force,&dist_ij,d1_theta);
+ v3_scale(&temp,&dist_ik,d_theta2);
+ v3_add(&force,&force,&temp);
+
+ /* part 1 of derivation of b_ij */
+ v3_scale(&force,sin_theta*2*h_cos*f_c_ik*frac1);
+
+ /* part 2 of derivation of b_ij */
+ v3_scale(&temp,&dist_ik,df_c_ik*bracket1);
+
+ /* sum up and scale ... */
+ v3_add(&temp,&temp,&force);
+ scale=bracket2_n_1*n*betan*(1+betan*bracket3_pow_1)*chi*(1.0/(2.0*n));
+ v3_scale(&temp,&temp,scale);
+
+ /* now construct an energy and a force out of that */
+ v3_scale(&temp,&temp,f_a);
+ v3_scale(&force,&dist_ij,df_a*b_ij);
+ v3_add(&temp,&temp,&force);
+ v3_scale(&temp,&temp,f_c);
+ v3_scale(&force,&dist_ij,df_c*b_ij*f_a);
+ v3_add(&force,&force,&temp);
+
+ /* add forces */
+ v3_add(&(ai->f),&(ai->f),&force);
+ /* energy is 0.5 f_r f_c, but we will sum it up twice ... */
+ moldyn->energy+=(0.25*f_a*b_ij*f_c);
+
return 0;
}
projects / physik / posic.git / blobdiff