#include "../moldyn.h"
#include "../math/math.h"
-//#include "tersoff.h"
+#include "tersoff.h"
/* create mixed terms from parameters and set them */
int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
t_tersoff_exchange *exchange;
brand=ai->brand;
- params=moldyn->pot1b_params;
+ params=moldyn->pot_params;
exchange=&(params->exchange);
/*
double s_r;
double arg;
- params=moldyn->pot2b_params;
+ /* use newtons third law */
+ //if(ai<aj) return 0;
+
+ params=moldyn->pot_params;
brand=aj->brand;
exchange=&(params->exchange);
*
*/
- /* constants */
- if(brand==ai->brand) {
- S=params->S[brand];
+ /* determine cutoff square */
+ if(brand==ai->brand)
S2=params->S2[brand];
- R=params->R[brand];
- A=params->A[brand];
- B=params->B[brand];
- lambda=params->lambda[brand];
- mu=params->mu[brand];
- exchange->chi=1.0;
- }
- else {
- S=params->Smixed;
+ else
S2=params->S2mixed;
- R=params->Rmixed;
- A=params->Amixed;
- B=params->Bmixed;
- lambda=params->lambda_m;
- mu=params->mu_m;
- params->exchange.chi=params->chi;
- }
/* dist_ij, d_ij */
v3_sub(&dist_ij,&(aj->r),&(ai->r));
exchange->d_j=&(params->d[brand]);
exchange->h_j=&(params->h[brand]);
if(brand==ai->brand) {
+ S=params->S[brand];
+ R=params->R[brand];
+ A=params->A[brand];
+ B=params->B[brand];
+ lambda=params->lambda[brand];
+ mu=params->mu[brand];
+ exchange->chi=1.0;
exchange->betajnj=exchange->betaini;
exchange->cj2=exchange->ci2;
exchange->dj2=exchange->di2;
exchange->cj2dj2=exchange->ci2di2;
}
else {
+ S=params->Smixed;
+ R=params->Rmixed;
+ A=params->Amixed;
+ B=params->Bmixed;
+ lambda=params->lambda_m;
+ mu=params->mu_m;
+ params->exchange.chi=params->chi;
exchange->betajnj=pow(*(exchange->beta_j),*(exchange->n_j));
exchange->cj2=params->c[brand]*params->c[brand];
exchange->dj2=params->d[brand]*params->d[brand];
v3_add(&(ai->f),&(ai->f),&force);
/* virial */
- ai->virial.xx-=force.x*dist_ij.x;
- ai->virial.yy-=force.y*dist_ij.y;
- ai->virial.zz-=force.z*dist_ij.z;
- ai->virial.xy-=force.x*dist_ij.y;
- ai->virial.xz-=force.x*dist_ij.z;
- ai->virial.yz-=force.y*dist_ij.z;
+ virial_calc(ai,&force,&dist_ij);
+ //ai->virial.xx-=force.x*dist_ij.x;
+ //ai->virial.yy-=force.y*dist_ij.y;
+ //ai->virial.zz-=force.z*dist_ij.z;
+ //ai->virial.xy-=force.x*dist_ij.y;
+ //ai->virial.xz-=force.x*dist_ij.z;
+ //ai->virial.yz-=force.y*dist_ij.z;
#ifdef DEBUG
if(ai==&(moldyn->atom[0])) {
double chi,ni,betaini,nj,betajnj;
double zeta;
- params=moldyn->pot2b_params;
+ params=moldyn->pot_params;
exchange=&(params->exchange);
/* we do not run if f_c_ij was detected to be 0! */
v3_add(&(ai->f),&(ai->f),&force);
/* virial */
- ai->virial.xx-=force.x*dist_ij->x;
- ai->virial.yy-=force.y*dist_ij->y;
- ai->virial.zz-=force.z*dist_ij->z;
- ai->virial.xy-=force.x*dist_ij->y;
- ai->virial.xz-=force.x*dist_ij->z;
- ai->virial.yz-=force.y*dist_ij->z;
+ virial_calc(ai,&force,dist_ij);
+ //ai->virial.xx-=force.x*dist_ij->x;
+ //ai->virial.yy-=force.y*dist_ij->y;
+ //ai->virial.zz-=force.z*dist_ij->z;
+ //ai->virial.xy-=force.x*dist_ij->y;
+ //ai->virial.xz-=force.x*dist_ij->z;
+ //ai->virial.yz-=force.y*dist_ij->z;
#ifdef DEBUG
if(ai==&(moldyn->atom[0])) {
/* virial - plus sign, as dist_ij = - dist_ji - (really??) */
// TEST ... with a minus instead
- ai->virial.xx-=force.x*dist_ij->x;
- ai->virial.yy-=force.y*dist_ij->y;
- ai->virial.zz-=force.z*dist_ij->z;
- ai->virial.xy-=force.x*dist_ij->y;
- ai->virial.xz-=force.x*dist_ij->z;
- ai->virial.yz-=force.y*dist_ij->z;
+ virial_calc(ai,&force,dist_ij);
+ //ai->virial.xx-=force.x*dist_ij->x;
+ //ai->virial.yy-=force.y*dist_ij->y;
+ //ai->virial.zz-=force.z*dist_ij->z;
+ //ai->virial.xy-=force.x*dist_ij->y;
+ //ai->virial.xz-=force.x*dist_ij->z;
+ //ai->virial.yz-=force.y*dist_ij->z;
#ifdef DEBUG
if(ai==&(moldyn->atom[0])) {
double tmp;
int brand;
- params=moldyn->pot3b_params;
+ params=moldyn->pot_params;
exchange=&(params->exchange);
if(!(exchange->run3bp))