// Hayter-Penfold (rescaled) MSA structure factor for screened Coulomb interactions // // C99 needs declarations of routines here double Iq(double QQ, double radius_effective, double zz, double VolFrac, double Temp, double csalt, double dialec); int sqcoef(int ir, double gMSAWave[]); int sqfun(int ix, int ir, double gMSAWave[]); double sqhcal(double qq, double gMSAWave[]); double Iq(double QQ, double radius_effective, double VolFrac, double zz, double Temp, double csalt, double dialec) { double gMSAWave[17]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17}; double Elcharge=1.602189e-19; // electron charge in Coulombs (C) double kB=1.380662e-23; // Boltzman constant in J/K double FrSpPerm=8.85418782E-12; //Permittivity of free space in C^2/(N m^2) double SofQ, Qdiam, Vp, ss; double SIdiam, diam, Kappa, cs, IonSt; double Perm, Beta; double charge; int ierr; diam=2*radius_effective; //in A //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////// convert to USEFUL inputs in SI units // //////////////////////////// NOTE: easiest to do EVERYTHING in SI units // //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// Beta=1.0/(kB*Temp); // in Joules^-1 Perm=dialec*FrSpPerm; //in C^2/(N m^2) charge=zz*Elcharge; //in Coulomb (C) SIdiam = diam*1.0E-10; //in m Vp=M_4PI_3*cube(SIdiam/2.0); //in m^3 cs=csalt*6.022E23*1.0E3; //# salt molecules/m^3 // Compute the derived values of : // Ionic strength IonSt (in C^2/m^3) // Kappa (Debye-Huckel screening length in m) // and gamma Exp(-k) // the zz*VolFrac/Vp is for the counterions from the micelle, assumed monovalent, the 2.0*cs if for added salt, assumed 1:1 electolyte IonSt=0.5 * Elcharge*Elcharge*(zz*VolFrac/Vp+2.0*cs); Kappa=sqrt(2*Beta*IonSt/Perm); //Kappa calc from Ionic strength // Kappa=2/SIdiam // Use to compare with HP paper gMSAWave[5]=Beta*charge*charge/(M_PI*Perm*SIdiam*square(2.0+Kappa*SIdiam)); // Finally set up dimensionless parameters Qdiam=QQ*diam; gMSAWave[6] = Kappa*SIdiam; gMSAWave[4] = VolFrac; //Function sqhpa(qq) {this is where Hayter-Penfold program began} // FIRST CALCULATE COUPLING ss=pow(gMSAWave[4],(1.0/3.0)); gMSAWave[9] = 2.0*ss*gMSAWave[5]*exp(gMSAWave[6]-gMSAWave[6]/ss); // CALCULATE COEFFICIENTS, CHECK ALL IS WELL // AND IF SO CALCULATE S(Q*SIG) ierr=0; ierr=sqcoef(ierr, gMSAWave); if (ierr>=0) { SofQ=sqhcal(Qdiam, gMSAWave); }else{ SofQ=NAN; // print "Error Level = ",ierr // print "Please report HPMSA problem with above error code" } return(SofQ); } ///////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////// // // // CALCULATES RESCALED VOLUME FRACTION AND CORRESPONDING // COEFFICIENTS FOR "SQHPA" // // JOHN B. HAYTER (I.L.L.) 14-SEP-81 // // ON EXIT: // // SETA IS THE RESCALED VOLUME FRACTION // SGEK IS THE RESCALED CONTACT POTENTIAL // SAK IS THE RESCALED SCREENING CONSTANT // A,B,C,F,U,V ARE THE MSA COEFFICIENTS // G1= G(1+) IS THE CONTACT VALUE OF G(R/SIG): // FOR THE GILLAN CONDITION, THE DIFFERENCE FROM // ZERO INDICATES THE COMPUTATIONAL ACCURACY. // // IR > 0: NORMAL EXIT, IR IS THE NUMBER OF ITERATIONS. // < 0: FAILED TO CONVERGE // int sqcoef(int ir, double gMSAWave[]) { int itm=40,ix,ig,ii; double acc=5.0E-6,del,e1,e2,f1,f2; // WAVE gMSAWave = $"root:HayPenMSA:gMSAWave" f1=0; //these were never properly initialized... f2=0; ig=1; if (gMSAWave[6]>=(1.0+8.0*gMSAWave[4])) { ig=0; gMSAWave[15]=gMSAWave[14]; gMSAWave[16]=gMSAWave[4]; ix=1; ir = sqfun(ix,ir,gMSAWave); gMSAWave[14]=gMSAWave[15]; gMSAWave[4]=gMSAWave[16]; if((ir<0.0) || (gMSAWave[14]>=0.0)) { return ir; } } gMSAWave[10]=fmin(gMSAWave[4],0.20); if ((ig!=1) || ( gMSAWave[9]>=0.15)) { ii=0; do { ii=ii+1; if(ii>itm) { ir=-1; return ir; } if (gMSAWave[10]<=0.0) { gMSAWave[10]=gMSAWave[4]/ii; } if(gMSAWave[10]>0.6) { gMSAWave[10] = 0.35/ii; } e1=gMSAWave[10]; gMSAWave[15]=f1; gMSAWave[16]=e1; ix=2; ir = sqfun(ix,ir,gMSAWave); f1=gMSAWave[15]; e1=gMSAWave[16]; e2=gMSAWave[10]*1.01; gMSAWave[15]=f2; gMSAWave[16]=e2; ix=2; ir = sqfun(ix,ir,gMSAWave); f2=gMSAWave[15]; e2=gMSAWave[16]; e2=e1-(e2-e1)*f1/(f2-f1); gMSAWave[10] = e2; del = fabs((e2-e1)/e1); } while (del>acc); gMSAWave[15]=gMSAWave[14]; gMSAWave[16]=e2; ix=4; ir = sqfun(ix,ir,gMSAWave); gMSAWave[14]=gMSAWave[15]; e2=gMSAWave[16]; ir=ii; if ((ig!=1) || (gMSAWave[10]>=gMSAWave[4])) { return ir; } } gMSAWave[15]=gMSAWave[14]; gMSAWave[16]=gMSAWave[4]; ix=3; ir = sqfun(ix,ir,gMSAWave); gMSAWave[14]=gMSAWave[15]; gMSAWave[4]=gMSAWave[16]; if ((ir>=0) && (gMSAWave[14]<0.0)) { ir=-3; } return ir; } int sqfun(int ix, int ir, double gMSAWave[]) { double acc=1.0e-6; double reta,eta2,eta21,eta22,eta3,eta32,eta2d,eta2d2,eta3d,eta6d,e12,e24,rgek; double rak,ak1,ak2,dak,dak2,dak4,d,d2,dd2,dd4,dd45,ex1,ex2,sk,ck,ckma,skma; double al1,al2,al3,al4,al5,al6,be1,be2,be3,vu1,vu2,vu3,vu4,vu5,ph1,ph2,ta1,ta2,ta3,ta4,ta5; double a1,a2,a3,b1,b2,b3,v1,v2,v3,p1,p2,p3,pp,pp1,pp2,p1p2,t1,t2,t3,um1,um2,um3,um4,um5,um6; double w0,w1,w2,w3,w4,w12,w13,w14,w15,w16,w24,w25,w26,w32,w34,w3425,w35,w3526,w36,w46,w56; double fa,fap,ca,e24g,pwk,qpw,pg,del,fun,fund,g24; int ii,ibig,itm=40; // WAVE gMSAWave = $"root:HayPenMSA:gMSAWave" a2=0; a3=0; b2=0; b3=0; v2=0; v3=0; p2=0; p3=0; // CALCULATE CONSTANTS; NOTATION IS HAYTER PENFOLD (1981) reta = gMSAWave[16]; eta2 = reta*reta; eta3 = eta2*reta; e12 = 12.0*reta; e24 = e12+e12; gMSAWave[13] = pow( (gMSAWave[4]/gMSAWave[16]),(1.0/3.0)); gMSAWave[12]=gMSAWave[6]/gMSAWave[13]; ibig=0; if (( gMSAWave[12]>15.0) && (ix==1)) { ibig=1; } gMSAWave[11] = gMSAWave[5]*gMSAWave[13]*exp(gMSAWave[6]- gMSAWave[12]); rgek = gMSAWave[11]; rak = gMSAWave[12]; ak2 = rak*rak; ak1 = 1.0+rak; dak2 = 1.0/ak2; dak4 = dak2*dak2; d = 1.0-reta; d2 = d*d; dak = d/rak; dd2 = 1.0/d2; dd4 = dd2*dd2; dd45 = dd4*2.0e-1; eta3d=3.0*reta; eta6d = eta3d+eta3d; eta32 = eta3+ eta3; eta2d = reta+2.0; eta2d2 = eta2d*eta2d; eta21 = 2.0*reta+1.0; eta22 = eta21*eta21; // ALPHA(I) al1 = -eta21*dak; al2 = (14.0*eta2-4.0*reta-1.0)*dak2; al3 = 36.0*eta2*dak4; // BETA(I) be1 = -(eta2+7.0*reta+1.0)*dak; be2 = 9.0*reta*(eta2+4.0*reta-2.0)*dak2; be3 = 12.0*reta*(2.0*eta2+8.0*reta-1.0)*dak4; // NU(I) vu1 = -(eta3+3.0*eta2+45.0*reta+5.0)*dak; vu2 = (eta32+3.0*eta2+42.0*reta-2.0e1)*dak2; vu3 = (eta32+3.0e1*reta-5.0)*dak4; vu4 = vu1+e24*rak*vu3; vu5 = eta6d*(vu2+4.0*vu3); // PHI(I) ph1 = eta6d/rak; ph2 = d-e12*dak2; // TAU(I) ta1 = (reta+5.0)/(5.0*rak); ta2 = eta2d*dak2; ta3 = -e12*rgek*(ta1+ta2); ta4 = eta3d*ak2*(ta1*ta1-ta2*ta2); ta5 = eta3d*(reta+8.0)*1.0e-1-2.0*eta22*dak2; // double PRECISION SINH(K), COSH(K) ex1 = exp(rak); ex2 = 0.0; if ( gMSAWave[12]<20.0) { ex2=exp(-rak); } sk=0.5*(ex1-ex2); ck = 0.5*(ex1+ex2); ckma = ck-1.0-rak*sk; skma = sk-rak*ck; // a(I) a1 = (e24*rgek*(al1+al2+ak1*al3)-eta22)*dd4; if (ibig==0) { a2 = e24*(al3*skma+al2*sk-al1*ck)*dd4; a3 = e24*(eta22*dak2-0.5*d2+al3*ckma-al1*sk+al2*ck)*dd4; } // b(I) b1 = (1.5*reta*eta2d2-e12*rgek*(be1+be2+ak1*be3))*dd4; if (ibig==0) { b2 = e12*(-be3*skma-be2*sk+be1*ck)*dd4; b3 = e12*(0.5*d2*eta2d-eta3d*eta2d2*dak2-be3*ckma+be1*sk-be2*ck)*dd4; } // V(I) v1 = (eta21*(eta2-2.0*reta+1.0e1)*2.5e-1-rgek*(vu4+vu5))*dd45; if (ibig==0) { v2 = (vu4*ck-vu5*sk)*dd45; v3 = ((eta3-6.0*eta2+5.0)*d-eta6d*(2.0*eta3-3.0*eta2+18.0*reta+1.0e1)*dak2+e24*vu3+vu4*sk-vu5*ck)*dd45; } // P(I) pp1 = ph1*ph1; pp2 = ph2*ph2; pp = pp1+pp2; p1p2 = ph1*ph2*2.0; p1 = (rgek*(pp1+pp2-p1p2)-0.5*eta2d)*dd2; if (ibig==0) { p2 = (pp*sk+p1p2*ck)*dd2; p3 = (pp*ck+p1p2*sk+pp1-pp2)*dd2; } // T(I) t1 = ta3+ta4*a1+ta5*b1; if (ibig!=0) { // VERY LARGE SCREENING: ASYMPTOTIC SOLUTION v3 = ((eta3-6.0*eta2+5.0)*d-eta6d*(2.0*eta3-3.0*eta2+18.0*reta+1.0e1)*dak2+e24*vu3)*dd45; t3 = ta4*a3+ta5*b3+e12*ta2 - 4.0e-1*reta*(reta+1.0e1)-1.0; p3 = (pp1-pp2)*dd2; b3 = e12*(0.5*d2*eta2d-eta3d*eta2d2*dak2+be3)*dd4; a3 = e24*(eta22*dak2-0.5*d2-al3)*dd4; um6 = t3*a3-e12*v3*v3; um5 = t1*a3+a1*t3-e24*v1*v3; um4 = t1*a1-e12*v1*v1; al6 = e12*p3*p3; al5 = e24*p1*p3-b3-b3-ak2; al4 = e12*p1*p1-b1-b1; w56 = um5*al6-al5*um6; w46 = um4*al6-al4*um6; fa = -w46/w56; ca = -fa; gMSAWave[3] = fa; gMSAWave[2] = ca; gMSAWave[1] = b1+b3*fa; gMSAWave[0] = a1+a3*fa; gMSAWave[8] = v1+v3*fa; gMSAWave[14] = -(p1+p3*fa); gMSAWave[15] = gMSAWave[14]; if (fabs(gMSAWave[15])<1.0e-3) { gMSAWave[15] = 0.0; } gMSAWave[10] = gMSAWave[16]; } else { t2 = ta4*a2+ta5*b2+e12*(ta1*ck-ta2*sk); t3 = ta4*a3+ta5*b3+e12*(ta1*sk-ta2*(ck-1.0))-4.0e-1*reta*(reta+1.0e1)-1.0; // MU(i) um1 = t2*a2-e12*v2*v2; um2 = t1*a2+t2*a1-e24*v1*v2; um3 = t2*a3+t3*a2-e24*v2*v3; um4 = t1*a1-e12*v1*v1; um5 = t1*a3+t3*a1-e24*v1*v3; um6 = t3*a3-e12*v3*v3; // GILLAN CONDITION ? // // YES - G(X=1+) = 0 // // COEFFICIENTS AND FUNCTION VALUE // if ((ix==1) || (ix==3)) { // NO - CALCULATE REMAINING COEFFICIENTS. // LAMBDA(I) al1 = e12*p2*p2; al2 = e24*p1*p2-b2-b2; al3 = e24*p2*p3; al4 = e12*p1*p1-b1-b1; al5 = e24*p1*p3-b3-b3-ak2; al6 = e12*p3*p3; // OMEGA(I) w16 = um1*al6-al1*um6; w15 = um1*al5-al1*um5; w14 = um1*al4-al1*um4; w13 = um1*al3-al1*um3; w12 = um1*al2-al1*um2; w26 = um2*al6-al2*um6; w25 = um2*al5-al2*um5; w24 = um2*al4-al2*um4; w36 = um3*al6-al3*um6; w35 = um3*al5-al3*um5; w34 = um3*al4-al3*um4; w32 = um3*al2-al3*um2; // w46 = um4*al6-al4*um6; w56 = um5*al6-al5*um6; w3526 = w35+w26; w3425 = w34+w25; // QUARTIC COEFFICIENTS w4 = w16*w16-w13*w36; w3 = 2.0*w16*w15-w13*w3526-w12*w36; w2 = w15*w15+2.0*w16*w14-w13*w3425-w12*w3526; w1 = 2.0*w15*w14-w13*w24-w12*w3425; w0 = w14*w14-w12*w24; // ESTIMATE THE STARTING VALUE OF f if (ix==1) { // LARGE K fap = (w14-w34-w46)/(w12-w15+w35-w26+w56-w32); } else { // ASSUME NOT TOO FAR FROM GILLAN CONDITION. // IF BOTH RGEK AND RAK ARE SMALL, USE P-W ESTIMATE. gMSAWave[14]=0.5*eta2d*dd2*exp(-rgek); if (( gMSAWave[11]<=2.0) && ( gMSAWave[11]>=0.0) && ( gMSAWave[12]<=1.0)) { e24g = e24*rgek*exp(rak); pwk = sqrt(e24g); qpw = (1.0-sqrt(1.0+2.0*d2*d*pwk/eta22))*eta21/d; gMSAWave[14] = -qpw*qpw/e24+0.5*eta2d*dd2; } pg = p1+gMSAWave[14]; ca = ak2*pg+2.0*(b3*pg-b1*p3)+e12*gMSAWave[14]*gMSAWave[14]*p3; ca = -ca/(ak2*p2+2.0*(b3*p2-b2*p3)); fap = -(pg+p2*ca)/p3; } // AND REFINE IT ACCORDING TO NEWTON ii=0; do { ii = ii+1; if (ii>itm) { // FAILED TO CONVERGE IN ITM ITERATIONS ir=-2; return (ir); } fa = fap; fun = w0+(w1+(w2+(w3+w4*fa)*fa)*fa)*fa; fund = w1+(2.0*w2+(3.0*w3+4.0*w4*fa)*fa)*fa; fap = fa-fun/fund; del=fabs((fap-fa)/fa); } while (del>acc); ir = ir+ii; fa = fap; ca = -(w16*fa*fa+w15*fa+w14)/(w13*fa+w12); gMSAWave[14] = -(p1+p2*ca+p3*fa); gMSAWave[15] = gMSAWave[14]; if (fabs(gMSAWave[15])<1.0e-3) { gMSAWave[15] = 0.0; } gMSAWave[10] = gMSAWave[16]; } else { ca = ak2*p1+2.0*(b3*p1-b1*p3); ca = -ca/(ak2*p2+2.0*(b3*p2-b2*p3)); fa = -(p1+p2*ca)/p3; if (ix==2) { gMSAWave[15] = um1*ca*ca+(um2+um3*fa)*ca+um4+um5*fa+um6*fa*fa; } if (ix==4) { gMSAWave[15] = -(p1+p2*ca+p3*fa); } } gMSAWave[3] = fa; gMSAWave[2] = ca; gMSAWave[1] = b1+b2*ca+b3*fa; gMSAWave[0] = a1+a2*ca+a3*fa; gMSAWave[8] = (v1+v2*ca+v3*fa)/gMSAWave[0]; } g24 = e24*rgek*ex1; gMSAWave[7] = (rak*ak2*ca-g24)/(ak2*g24); return (ir); } double sqhcal(double qq, double gMSAWave[]) { double SofQ,etaz,akz,gekz,e24,x1,x2,ck,sk,ak2,qk,q2k,qk2,qk3,qqk,sink,cosk,asink,qcosk,aqk,inter; // WAVE gMSAWave = $"root:HayPenMSA:gMSAWave" etaz = gMSAWave[10]; akz = gMSAWave[12]; gekz = gMSAWave[11]; e24 = 24.0*etaz; x1 = exp(akz); x2 = 0.0; if ( gMSAWave[12]<20.0) { x2 = exp(-akz); } ck = 0.5*(x1+x2); sk = 0.5*(x1-x2); ak2 = akz*akz; qk = qq/gMSAWave[13]; q2k = qk*qk; if (qk<=1.0e-08) { SofQ = -1.0/gMSAWave[0]; } else { // this rescales Q.sigma = 2.Q.Radius, so is hard to predict the value to test the function if (qk<=0.01) { // try Taylor series expansion at small qk (RKH Feb 2016, with help from Mathematica), // transition point may need to depend on precision of cpu used and ALAS on the values of some of the parameters ! // note have adsorbed a factor 24 from SofQ= // needs thorough test over wide range of parameter space! // there seem to be some rounding issues here in single precision, must use double aqk = gMSAWave[0]*(8.0+2.0*etaz) + 6*gMSAWave[1] -12.0*gMSAWave[3] -24*(gekz*(1.0+akz) -ck*akz*gMSAWave[2] +gMSAWave[3]*(ck-1.0) +(gMSAWave[2]-gMSAWave[3]*akz)*sk )/ak2 +q2k*( -(gMSAWave[0]*(48.0+15.0*etaz) +40.0*gMSAWave[1])/60.0 +gMSAWave[3] +(4.0/ak2)*(gekz*(9.0+7.0*akz) +ck*(9.0*gMSAWave[3] -7.0*gMSAWave[2]*akz) +sk*(9.0*gMSAWave[2] -7.0*gMSAWave[3]*akz)) ); SofQ = 1.0/(1.0-gMSAWave[10]*aqk); } else { qk2 = 1.0/q2k; qk3 = qk2/qk; qqk = 1.0/(qk*(q2k+ak2)); SINCOS(qk,sink,cosk); asink = akz*sink; qcosk = qk*cosk; aqk = gMSAWave[0]*(sink-qcosk); aqk=aqk+gMSAWave[1]*((2.0*qk2-1.0)*qcosk+2.0*sink-2.0/qk); inter=24.0*qk3+4.0*(1.0-6.0*qk2)*sink; aqk=(aqk+0.5*etaz*gMSAWave[0]*(inter-(1.0-12.0*qk2+24.0*qk2*qk2)*qcosk))*qk3; aqk=aqk +gMSAWave[2]*(ck*asink-sk*qcosk)*qqk; aqk=aqk +gMSAWave[3]*(sk*asink-qk*(ck*cosk-1.0))*qqk; aqk=aqk +gMSAWave[3]*(cosk-1.0)*qk2; aqk=aqk -gekz*(asink+qcosk)*qqk; SofQ = 1.0/(1.0 -e24*aqk); } } return (SofQ); }