C Computer program to calculate the reorganization free energy for electron transfer at the liquid/liquid interface. An analytical solution for the one-dimensional case

#include stdio.h

#include math.h

#define PI 3.14159265358979323846

#define abs(x) ((x)>0.?(x):-(x))

#define sgn(x) ((x)>0.?1.:-1.)

#define NZ 80

double ra,rb;

main(argc,argv)

int argc;

char *argv[];

{ double Ir();

double z,Rab,e1s,e2s,e1o,e2o,Ha,Hb,ES,ES0;

double rA, rB;

int i,j,loc;

char buf[80];

FILE *fp;

if (argc != 8){

fprintf(stderr,"Usage: %s ra rb Rab e1s e2s e1o e2o \n", argv[0]);

exit(1);

} rA = atof(argv[1]);

rB = atof(argv[2]);

Rab = atof(argv[3]);

e1s = atof(argv[4]);

e2s = atof(argv[5]);

e1o = atof(argv[6]);

e2o = atof(argv[7]);

if (Rab < rA+rB){

fprintf(stderr,"Rab = %f must be greater than rA+rB = %f\n",Rab,rA+rB);

exit(1);

} ES0 = 0.5*(1.0/rA+1.0/rB)*(1.0/e1o-1.0/e1s);

/*Main loop over the location of the center point of the ion-pair*/ printf("z/R ");

for (j=1;

j<41;j++){

Rab = 4+j*0.05;

printf("%-7.2f",Rab);

} printf("\n");

for (i=0;

i z = -8+16.0*(i+0.333)/(1.0*NZ);

printf("%-7.2f",z);

for (j=0;

j<40;j++){

Rab = 4+j*0.05;

Ha = abs(z-0.5*Rab);

Hb = abs(z+0.5*Rab);

if (z <= -0.5*Rab-rB) loc = 1;

if (z <-0.5*Rab && z > -0.5*Rab-rB) loc = 2;

if (z <= -0.5*Rab+rB && z > -0.5*Rab) loc = 3;

if (z <= 0.5*Rab-rA && z > -0.5*Rab+rB) loc = 4;

if (z <= 0.5*Rab && z > 0.5*Rab-rA) loc = 5;

if (z <= 0.5*Rab+rA && z > 0.5*Rab) loc = 6;

if (z > 0.5*Rab+rA ) loc = 7;

if (loc <5){

ra = rA;

rb = rB;

ES = Ir(e1o,e2o,Ha,Hb,loc)-Ir(e1s,e2s,Ha,Hb,loc);

} if (loc > 4){

ra = rB;

rb = rA;

ES = Ir(e2o,e1o,Hb,Ha,8-loc)-Ir(e2s,e1s,Hb,Ha,8-loc);

} printf("%-7.3f",ES/ES0);

} printf("\n");

} } double Ir(e1,e2,ha,hb,n)

double e1,e2,ha,hb;

int n;

{ double eta,val,valA,valB,valB1,valB2,sint();

eta = (e1-e2)/(e1+e2);

switch (n){

case 1:

/*both ions are totaly in the e1 region (z<0)*/

valA = sint(-ra,ra,0.,0.,ra) + eta*sint(-ra,ra,2*ha,0.,ra) - sint(-ra,ra,ha-hb,0.,ra) - eta*sint(-ra,ra,ha+hb,0.,ra) + eta*sint(-ra,ra,0.,2*ha,ra) + eta*eta*sint(-ra,ra,2*ha,2*ha,ra) - eta*sint(-ra,ra,ha-hb,2*ha,ra) - eta*eta*sint(-ra,ra,ha+hb,2*ha,ra) - sint(-ra,ra,0.,ha-hb,ra) - eta*sint(-ra,ra,2*ha,ha-hb,ra) + sint(-ra,ra,ha-hb,ha-hb,ra) + eta*sint(-ra,ra,ha+hb,ha-hb,ra) - eta*sint(-ra,ra,0.,ha+hb,ra) - eta*eta*sint(-ra,ra,2*ha,ha+hb,ra) + eta*sint(-ra,ra,ha-hb,ha+hb,ra) + eta*eta*sint(-ra,ra,ha+hb,ha+hb,ra);

valA /= 8*PI*e1;

valB = sint(-rb,rb,hb-ha,hb-ha,rb) + eta*sint(-rb,rb,ha+hb,hb-ha,rb) - sint(-rb,rb,0.,hb-ha,rb) - eta*sint(-rb,rb,2*hb,hb-ha,rb) + eta*sint(-rb,rb,hb-ha,ha+hb,rb) + eta*eta*sint(-rb,rb,ha+hb,ha+hb,rb) - eta*sint(-rb,rb,0.,ha+hb,rb) - eta*eta*sint(-rb,rb,2*hb,ha+hb,rb) - sint(-rb,rb,hb-ha,0.,rb) - eta*sint(-rb,rb,ha+hb,0.,rb) + sint(-rb,rb,0.,0.,rb) + eta*sint(-rb,rb,2*hb,0.,rb) - eta*sint(-rb,rb,hb-ha,2*hb,rb) - eta*eta*sint(-rb,rb,ha+hb,2*hb,rb) + eta*sint(-rb,rb,0.,2*hb,rb) + eta*eta*sint(-rb,rb,2*hb,2*hb,rb);

valB /= 8*PI*e1;

val = valA + valB;

break;

case 2:

/*both centers are in e1, one ion is in e1 2nd ion is partially in e1*/

valA /= 8*PI*e1;

valB1 = sint(-hb,rb,hb-ha,hb-ha,rb) + eta*sint(-hb,rb,ha+hb,hb-ha,rb) - sint(-hb,rb,0.,hb-ha,rb) - eta*sint(-hb,rb,2*hb,hb-ha,rb) + eta*sint(-hb,rb,hb-ha,ha+hb,rb) + eta*eta*sint(-hb,rb,ha+hb,ha+hb,rb) - eta*sint(-hb,rb,0.,ha+hb,rb) - eta*eta*sint(-hb,rb,2*hb,ha+hb,rb) - sint(-hb,rb,hb-ha,0.,rb) - eta*sint(-hb,rb,ha+hb,0.,rb) + sint(-hb,rb,0.,0.,rb) + eta*sint(-hb,rb,2*hb,0.,rb) - eta*sint(-hb,rb,hb-ha,2*hb,rb) - eta*eta*sint(-hb,rb,ha+hb,2*hb,rb) + eta*sint(-hb,rb,0.,2*hb,rb) + eta*eta*sint(-hb,rb,2*hb,2*hb,rb);

valB1 /= 8*PI*e1;

valB2 = sint(-rb,-hb,hb-ha,hb-ha,rb) - sint(-rb,-hb,hb-ha,0.,rb) - sint(-rb,-hb,0.,hb-ha,rb) + sint(-rb,-hb,0.,0.,rb);

valB2 *= (e2/(2*PI*(e1+e2)*(e1+e2)));

val = valA + valB1 + valB2;

break;

case 3:

/*one ion is in e1 2nd ion is partially in e1, but its center is in e2*/

valA = sint(-ra,ra,0.,0.,ra) + eta*sint(-ra,ra,2*ha,0.,ra) - (1+eta)*sint(-ra,ra,ha+hb,0.,ra) + eta*sint(-ra,ra,0.,2*ha,ra) + eta*eta*sint(-ra,ra,2*ha,2*ha,ra) - eta*(1+eta)*sint(-ra,ra,ha+hb,2*ha,ra) - (1+eta)*sint(-ra,ra,0.,ha+hb,ra) - eta*(1+eta)*sint(-ra,ra,2*ha,ha+hb,ra) + (1+eta)*(1+eta)*sint(-ra,ra,ha+hb,ha+hb,ra);

valA /= 8*PI*e1;

valB1 = sint(hb,rb,-ha-hb,-ha-hb,rb) + eta*sint(hb,rb,ha-hb,-ha-hb,rb) - (1+eta)*sint(hb,rb,0.,-ha-hb,rb) + eta*sint(hb,rb,-ha-hb,ha-hb,rb) + eta*eta*sint(hb,rb,ha-hb,ha-hb,rb) - eta*(1+eta)*sint(hb,rb,0.,ha-hb,rb) - (1+eta)*sint(hb,rb,-ha-hb,0.,rb) - eta*(1+eta)*sint(hb,rb,ha-hb,0.,rb) + (1+eta)*(1+eta)*sint(hb,rb,0.,0.,rb);

valB1 /= 8*PI*e1;

valB2 = sint(-rb,hb,0.,0.,rb) - eta*sint(-rb,hb,-2*hb,0.,rb) - (1-eta)*sint(-rb,hb,-ha-hb,0.,rb) - eta*sint(-rb,hb,0.,-2*hb,rb) + eta*eta*sint(-rb,hb,-2*hb,-2*hb,rb) + eta*(1-eta)*sint(-rb,hb,-ha-hb,-2*hb,rb) - (1-eta)*sint(-rb,hb,0.,-ha-hb,rb) + eta*(1-eta)*sint(-rb,hb,-2*hb,-ha-hb,rb) + (1-eta)*(1-eta)*sint(-rb,hb,-ha-hb,-ha-hb,rb);

valB2 /= 8*PI*e2;

val = valA + valB1 + valB2;

break;

case 4:

/*one ion is totaly in e1 2nd ion is totaly in e2*/

valA /= 8*PI*e1;

valB = sint(-rb,rb,0.,0.,rb) - eta*sint(-rb,rb,-2*hb,0.,rb) - (1-eta)*sint(-rb,rb,-ha-hb,0.,rb) - eta*sint(-rb,rb,0.,-2*hb,rb) + eta*eta*sint(-rb,rb,-2*hb,-2*hb,rb) + eta*(1-eta)*sint(-rb,rb,-ha-hb,-2*hb,rb) - (1-eta)*sint(-rb,rb,0.,-ha-hb,rb) + eta*(1-eta)*sint(-rb,rb,-2*hb,-ha-hb,rb) + (1-eta)*(1-eta)*sint(-rb,rb,-ha-hb,-ha-hb,rb);

valB /= 8*PI*e2;

val = valA + valB;

break;

} return(val);

} double sint(zmin,zmax,alpha,beta,r)

double zmin,zmax,alpha,beta,r;

{ double ab,A,B,Ap,Bp,Am,Bm,val;

ab = alpha*beta;

if (ab != 0){

A = (r*r+alpha*alpha)/(2*abs(alpha));

B = (r*r+beta*beta)/(2*abs(beta));

Ap = sqrt(A+sgn(alpha)*zmax);

Bp = sqrt(B+sgn(beta)*zmax);

Am = sqrt(A+sgn(alpha)*zmin);

Bm = sqrt(B+sgn(beta)*zmin);

if (ab >0){/* both positive or both negative*/

val = log((Ap+Bp)/(Am+Bm)) -((r*r-beta*beta)/(2*abs(beta)))*(1.0/(Bp*(Ap+Bp))- 1.0/(Bm*(Am+Bm)));

val *= sgn(alpha)*PI/sqrt(ab);

} if (ab < 0){

val = - (r*r-beta*beta)/(2*beta*(A+B))*(Ap/Bp-Am/Bm);

if (alpha > 0){

val += atan(Ap/Bp)-atan(Am/Bm);

val *= PI/sqrt(-ab);

} if (alpha < 0){

val += atan(Bp/Ap)-atan(Bm/Am);

val *= PI/sqrt(-ab);

} } } else {/*one or both alpha and beta are zero*/

if (beta == 0 && alpha != 0){

val = sqrt(r*r+alpha*alpha+2*alpha*zmax)-sqrt(r*r+alpha*alpha+2*alpha*zmin);

val *= (2*PI/(alpha*r));

} if (alpha == 0 && beta != 0){

val = (beta+zmax)/sqrt(r*r+beta*beta+2*beta*zmax);

val -= (beta+zmin)/sqrt(r*r+beta*beta+2*beta*zmin);

val *= (2*PI/r);

} if (alpha == 0 && beta == 0) val = 2*PI*(zmax-zmin)/(r*r);

} return(val);