--- /dev/null
+/*
+ * albe.c - albe potential
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <math.h>
+
+#include "../moldyn.h"
+#include "../math/math.h"
+#include "albe.h"
+
+/* create mixed terms from parameters and set them */
+int albe_mult_complete_params(t_albe_mult_params *p) {
+
+ printf("[moldyn] albe parameter completion\n");
+ p->S2[0]=p->S[0]*p->S[0];
+ p->S2[1]=p->S[1]*p->S[1];
+ p->S2mixed=p->Smixed*p->Smixed;
+
+ printf("[moldyn] albe mult parameter info:\n");
+ printf(" S (A) | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
+ printf(" R (A) | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
+ printf(" A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
+ printf(" B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
+ printf(" lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
+ p->lambda_m);
+ printf(" mu | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
+ printf(" gamma | %f | %f\n",p->gamma[0],p->gamma[1]);
+ printf(" c | %f | %f\n",p->c[0],p->c[1]);
+ printf(" d | %f | %f\n",p->d[0],p->d[1]);
+ printf(" h | %f | %f\n",p->h[0],p->h[1]);
+
+ return 0;
+}
+
+/* albe 1 body part */
+int albe_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
+
+ int brand;
+ t_albe_mult_params *params;
+ t_albe_exchange *exchange;
+
+ brand=ai->brand;
+ params=moldyn->pot_params;
+ exchange=&(params->exchange);
+
+ /*
+ * simple: point constant parameters only depending on atom i to
+ * their right values
+ */
+
+ exchange->gamma_i=&(params->gamma[brand]);
+ exchange->c_i=&(params->c[brand]);
+ exchange->d_i=&(params->d[brand]);
+ exchange->h_i=&(params->h[brand]);
+
+ exchange->ci2=params->c[brand]*params->c[brand];
+ exchange->di2=params->d[brand]*params->d[brand];
+ exchange->ci2di2=exchange->ci2/exchange->di2;
+
+ return 0;
+}
+
+/* albe 3 body potential function (first ij loop) */
+int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+
+ t_albe_mult_params *params;
+ t_albe_exchange *exchange;
+ unsigned char brand;
+ double S2;
+ t_3dvec dist_ij;
+ double d_ij2,d_ij;
+
+ params=moldyn->pot_params;
+ exchange=&(params->exchange);
+
+ /* reset zeta sum */
+ exchange->zeta_ij=0.0;
+
+ /*
+ * set ij depending values
+ */
+
+ brand=ai->brand;
+
+ if(brand==aj->brand)
+ S2=params->S2[brand];
+ else
+ S2=params->S2mixed;
+
+ /* dist_ij, d_ij2 */
+ v3_sub(&dist_ij,&(aj->r),&(ai->r));
+ if(bc) check_per_bound(moldyn,&dist_ij);
+ d_ij2=v3_absolute_square(&dist_ij);
+
+ /* if d_ij2 > S2 => no force & potential energy contribution */
+ if(d_ij2>S2) {
+ moldyn->run3bp=0;
+ return 0;
+ }
+
+ /* d_ij */
+ d_ij=sqrt(d_ij2);
+
+ /* store values */
+ exchange->dist_ij=dist_ij;
+ exchange->d_ij2=d_ij2;
+ exchange->d_ij=d_ij;
+
+ /* reset k counter for first k loop */
+ exchange->kcount=0;
+
+ return 0;
+}
+
+/* albe 3 body potential function (first k loop) */
+int albe_mult_3bp_k1(t_moldyn *moldyn,
+ t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+
+ t_albe_mult_params *params;
+ t_albe_exchange *exchange;
+ unsigned char brand;
+ double R,S,S2;
+ t_3dvec dist_ij,dist_ik;
+ double d_ik2,d_ik,d_ij;
+ double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
+ double f_c_ik,df_c_ik;
+ int kcount;
+
+ params=moldyn->pot_params;
+ exchange=&(params->exchange);
+ kcount=exchange->kcount;
+
+ if(kcount>ALBE_MAXN) {
+ printf("FATAL: neighbours = %d\n",kcount);
+ printf(" -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
+ }
+
+ /* ik constants */
+ brand=ai->brand;
+ if(brand==ak->brand) {
+ R=params->R[brand];
+ S=params->S[brand];
+ S2=params->S2[brand];
+ }
+ else {
+ R=params->Rmixed;
+ S=params->Smixed;
+ S2=params->S2mixed;
+ }
+
+ /* dist_ik, d_ik2 */
+ v3_sub(&dist_ik,&(ak->r),&(ai->r));
+ if(bc) check_per_bound(moldyn,&dist_ik);
+ d_ik2=v3_absolute_square(&dist_ik);
+
+ /* store data for second k loop */
+ exchange->dist_ik[kcount]=dist_ik;
+ exchange->d_ik2[kcount]=d_ik2;
+
+ /* return if not within cutoff */
+ if(d_ik2>S2) {
+ exchange->kcount++;
+ return 0;
+ }
+
+ /* d_ik */
+ d_ik=sqrt(d_ik2);
+
+ /* dist_ij, d_ij */
+ dist_ij=exchange->dist_ij;
+ d_ij=exchange->d_ij;
+
+ /* cos theta */
+ cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
+
+ /* g_ijk */
+ h_cos=*(exchange->h_i)-cos_theta;
+ d2_h_cos2=exchange->di2+(h_cos*h_cos);
+ frac=exchange->ci2/d2_h_cos2;
+ g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
+ dg=-2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2;
+
+ /* zeta sum += f_c_ik * g_ijk */
+ if(d_ik<=R) {
+ exchange->zeta_ij+=g;
+ f_c_ik=1.0;
+ df_c_ik=0.0;
+ }
+ else {
+ s_r=S-R;
+ arg=M_PI*(d_ik-R)/s_r;
+ f_c_ik=0.5+0.5*cos(arg);
+ df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
+ exchange->zeta_ij+=f_c_ik*g;
+ }
+
+ /* store even more data for second k loop */
+ exchange->g[kcount]=g;
+ exchange->dg[kcount]=dg;
+ exchange->d_ik[kcount]=d_ik;
+ exchange->cos_theta[kcount]=cos_theta;
+ exchange->f_c_ik[kcount]=f_c_ik;
+ exchange->df_c_ik[kcount]=df_c_ik;
+
+ /* increase k counter */
+ exchange->kcount++;
+
+ return 0;
+}
+
+int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+
+ t_albe_mult_params *params;
+ t_albe_exchange *exchange;
+ t_3dvec force;
+ double f_a,df_a,b,db,f_c,df_c;
+ double f_r,df_r;
+ double scale;
+ double mu,B;
+ double lambda,A;
+ double d_ij,r0;
+ unsigned char brand;
+ double S,R,s_r,arg;
+
+ params=moldyn->pot_params;
+ exchange=&(params->exchange);
+
+ brand=aj->brand;
+ if(brand==ai->brand) {
+ S=params->S[brand];
+ R=params->R[brand];
+ B=params->B[brand];
+ A=params->A[brand];
+ r0=params->r0[brand];
+ mu=params->mu[brand];
+ lambda=params->lambda[brand];
+ }
+ else {
+ S=params->Smixed;
+ R=params->Rmixed;
+ B=params->Bmixed;
+ A=params->Amixed;
+ r0=params->r0_mixed;
+ mu=params->mu_m;
+ lambda=params->lambda_m;
+ }
+
+ d_ij=exchange->d_ij;
+
+ /* f_c, df_c */
+ if(d_ij<R) {
+ f_c=1.0;
+ df_c=0.0;
+ }
+ else {
+ s_r=S-R;
+ arg=M_PI*(d_ij-R)/s_r;
+ f_c=0.5+0.5*cos(arg);
+ df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
+ }
+
+ /* f_a, df_a */
+ f_a=-B*exp(-mu*(d_ij-r0));
+ df_a=mu*f_a/d_ij;
+
+ /* f_r, df_r */
+ f_r=A*exp(-lambda*(d_ij-r0));
+ df_r=lambda*f_r/d_ij;
+
+ /* b, db */
+ if(exchange->zeta_ij==0.0) {
+ b=1.0;
+ db=0.0;
+ }
+ else {
+ b=1.0/sqrt(1.0+exchange->zeta_ij);
+ db=-0.5*b/(1.0+exchange->zeta_ij);
+ }
+
+ /* force contribution */
+ scale=-0.5*(f_c*(df_r+b*df_a)+df_c*(f_r+b*df_a));
+ v3_scale(&force,&(exchange->dist_ij),scale);
+ v3_add(&(ai->f),&(ai->f),&force);
+ v3_sub(&(aj->f),&(aj->f),&force); // dri rij = - drj rij
+
+#ifdef DEBUG
+ if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) {
+ printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
+ printf("adding %f %f %f\n",force.x,force.y,force.z);
+ if(ai==&(moldyn->atom[0]))
+ printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
+ if(aj==&(moldyn->atom[0]))
+ printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
+ printf("energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
+ f_c,b,f_a,f_r);
+ }
+#endif
+
+ /* virial */
+ if(aj<ai)
+ virial_calc(ai,&force,&(exchange->dist_ij));
+
+ /* dzeta prefactor = - 0.5 f_c f_a db */
+ exchange->pre_dzeta=-0.5*f_a*f_c*db;
+
+ /* energy contribution */
+ moldyn->energy+=0.5*f_c*(f_r+b*f_a);
+
+ /* reset k counter for second k loop */
+ exchange->kcount=0;
+
+ return 0;
+}
+
+/* albe 3 body potential function (second k loop) */
+int albe_mult_3bp_k2(t_moldyn *moldyn,
+ t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+
+ t_albe_mult_params *params;
+ t_albe_exchange *exchange;
+ int kcount;
+ t_3dvec dist_ik,dist_ij;
+ double d_ik2,d_ik,d_ij2,d_ij;
+ unsigned char brand;
+ double S2;
+ double g,dg,cos_theta;
+ double pre_dzeta;
+ double f_c_ik,df_c_ik;
+ double dijdik_inv,fcdg,dfcg;
+ t_3dvec dcosdri,dcosdrj,dcosdrk;
+ t_3dvec force,tmp;
+
+ params=moldyn->pot_params;
+ exchange=&(params->exchange);
+ kcount=exchange->kcount;
+
+ if(kcount>ALBE_MAXN)
+ printf("FATAL: neighbours!\n");
+
+ /* d_ik2 */
+ d_ik2=exchange->d_ik2[kcount];
+
+ brand=ak->brand;
+ if(brand==ai->brand)
+ S2=params->S2[brand];
+ else
+ S2=params->S2mixed;
+
+ /* return if d_ik > S */
+ if(d_ik2>S2) {
+ exchange->kcount++;
+ return 0;
+ }
+
+ /* prefactor dzeta */
+ pre_dzeta=exchange->pre_dzeta;
+
+ /* dist_ik, d_ik */
+ dist_ik=exchange->dist_ik[kcount];
+ d_ik=exchange->d_ik[kcount];
+
+ /* f_c_ik, df_c_ik */
+ f_c_ik=exchange->f_c_ik[kcount];
+ df_c_ik=exchange->df_c_ik[kcount];
+
+ /* dist_ij, d_ij2, d_ij */
+ dist_ij=exchange->dist_ij;
+ d_ij2=exchange->d_ij2;
+ d_ij=exchange->d_ij;
+
+ /* g, dg, cos_theta */
+ g=exchange->g[kcount];
+ dg=exchange->dg[kcount];
+ cos_theta=exchange->cos_theta[kcount];
+
+ /* cos_theta derivatives wrt i,j,k */
+ dijdik_inv=1.0/(d_ij*d_ik);
+ v3_scale(&dcosdrj,&dist_ik,dijdik_inv);
+ v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
+ v3_add(&dcosdrj,&dcosdrj,&tmp);
+ v3_scale(&dcosdrk,&dist_ij,dijdik_inv);
+ v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
+ v3_add(&dcosdrk,&dcosdrk,&tmp);
+ v3_add(&dcosdri,&dcosdrj,&dcosdrk);
+ v3_scale(&dcosdri,&dcosdri,-1.0);
+
+ /* f_c_ik * dg, df_c_ik * g */
+ fcdg=f_c_ik*dg;
+ dfcg=df_c_ik*g;
+
+ /* derivative wrt i */
+ v3_scale(&force,&dist_ik,dfcg);
+ v3_scale(&tmp,&dcosdri,fcdg);
+ v3_add(&force,&force,&tmp);
+ v3_scale(&force,&force,pre_dzeta);
+
+ /* force contribution */
+ v3_add(&(ai->f),&(ai->f),&force);
+
+#ifdef DEBUG
+ if(ai==&(moldyn->atom[0])) {
+ printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
+ printf("adding %f %f %f\n",force.x,force.y,force.z);
+ printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
+ }
+#endif
+
+ /* virial */
+ //virial_calc(ai,&force,&dist_ij);
+
+ /* derivative wrt j */
+ v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
+
+ /* force contribution */
+ v3_add(&(aj->f),&(aj->f),&force);
+
+#ifdef DEBUG
+ if(aj==&(moldyn->atom[0])) {
+ printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
+ printf("adding %f %f %f\n",force.x,force.y,force.z);
+ printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
+ }
+#endif
+
+ /* virial */
+ //v3_scale(&force,&force,-1.0);
+ if(aj<ai)
+ virial_calc(ai,&force,&dist_ij);
+
+ /* derivative wrt k */
+ v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
+ v3_scale(&tmp,&dcosdrk,fcdg);
+ v3_add(&force,&force,&tmp);
+ v3_scale(&force,&force,pre_dzeta);
+
+ /* force contribution */
+ v3_add(&(ak->f),&(ak->f),&force);
+
+#ifdef DEBUG
+ if(ak==&(moldyn->atom[0])) {
+ printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
+ printf("adding %f %f %f\n",force.x,force.y,force.z);
+ printf("total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
+ }
+#endif
+
+ /* virial */
+ //v3_scale(&force,&force,-1.0);
+ if(aj<ai)
+ virial_calc(ai,&force,&dist_ik);
+
+ /* increase k counter */
+ exchange->kcount++;
+
+ return 0;
+
+}
projects / physik / posic.git / blobdiff