PROGRAM PESDIATATDIAT IMPLICIT REAL*8 (A-H,O-Z) parameter(pigreco = 3.1415926535897932d0) COMMON/PARAM2/epsi1,epsi2,rm1,rm2 OPEN(7,file='prova.dat') r=2.2d0 do 10 i=1,700 c------------------------------------------------------- c He c------------------------------------------- c CS rrh2=1.5346d0 c #*1.15d0 c------------------------------------------- thetaa=90.d0 thetarada = thetaa*pigreco/180.0d0 CALL SUPATATNEW(R,rrh2,THETARADA,V1) pot90=V1 c--------------------------------------------- thetaa=0.d0 thetarada = thetaa*pigreco/180.0d0 CALL SUPATATNEW(R,rrh2,THETARADA,V1) pot0=V1 c--------------------------------------------- thetaa=180.d0 thetarada = thetaa*pigreco/180.0d0 CALL SUPATATNEW(R,rrh2,THETARADA,V1) pot180=V1 c--------------------------------------------- thetaa=30.d0 thetarada = thetaa*pigreco/180.0d0 CALL SUPATATNEW(R,rrh2,THETARADA,V1) pot30=V1 c--------------------------------------------- thetaa=60.d0 thetarada = thetaa*pigreco/180.0d0 CALL SUPATATNEW(R,rrh2,THETARADA,V1) pot60=V1 c--------------------------------------------- thetaa=120.d0 thetarada = thetaa*pigreco/180.0d0 CALL SUPATATNEW(R,rrh2,THETARADA,V1) pot120=V1 c--------------------------------------------- thetaa=150.d0 thetarada = thetaa*pigreco/180.0d0 CALL SUPATATNEW(R,rrh2,THETARADA,V1) pot150=V1 write(7,3) r,pot0,pot30,pot60,pot90,pot120,pot150,pot180 write(*,3) r,pot0,pot30,pot60,pot90,pot120,pot150,pot180 r=r+0.01d0 10 continue 3 format(8(1x,f11.4)) 4 format(f5.2,36(1x,f11.8)) stop end c-------------------------------------------------------------------------- c FROM HERE IF YOU USE ANOTHER PROGRAM ---------------- c-------------------------------------------------------------------------- SUBROUTINE SUPATATNEW(R,rr2,THETARADB,VTOT) IMPLICIT REAL*8 (A-H,O-Z) c R distanza intermolecolare in Angstr. (distanza tra He e centro di massa CS) c rr2 distanza interna (r_CS) molecola CS in angstr c THETARADB angolo di jacobi tra asse diatomo e distanza intermolecolare in radianti c c valore distanza equilibrio r_CS=1.5346 angstr. c VTOT in meV COMMON/PARAM2/epsi1,epsi2,rm1,rm2 COMMON/PARAM5/dist dimension xxx(6),yyy(6),zzz(6),dist(42),q(6) integer hh parameter(pigreco = 3.1415926535897932d0) real*8 nlj,mlj c prima specie (He) c polarizzabilita' (atomic units) alfaA=1.35d0 c------------------------------------------------- c nuova parte elettrostatica per lo stretching c------------------------------------------------- c seconda molecola (CS) c calcolo MAXBART aug-cc-pV5Z spazio occ,8,4,4,0; ACPF c Quadrupolo QB=-1.749141d0+2.5d0*(rr2-1.5346d0) #-1.d0*(rr2-1.5346d0)**2 #+0.25d0*(rr2-1.5346d0)**3 #-1.2d0*(rr2-1.5346d0)**4 c cariche parziali qbpos=QB/(2.d0*(rr2/0.529177d0/2.d0)**2) c write(*,*) QB,qbpos c stop c valore relativo a distribuzione CILINDRICA qbneg=-qbpos/2.d0 c Dipolo DIPB=0.775837d0-1.45d0*(rr2-1.5346d0) #+0.3d0*(rr2-1.5346d0)**2 #+0.45d0*(rr2-1.5346d0)**3 #-0.17*(rr2-1.5346d0)**4 c distanza equilibrio CS rrcs=1.5346d0 c coordinate del CM dell'He xxx(3)=R*dsin(thetaradb) yyy(3)=0.d0 zzz(3)=R*dcos(thetaradb) c fissa le coordinate del diatomo 2 rrC=rr2*32./44. rrS=rr2-rrC xxx(4)=0.d0 xxx(5)=0.d0 yyy(4)=0.d0 yyy(5)=0.d0 zzz(4)=rrS ! S zzz(5)=-rrC ! C c coordinate del CM di CS xxx(6)=0.d0 yyy(6)=0.d0 zzz(6)=0.d0 c cariche su CS c--------------------------------------------------------------------- q(4)=qbpos q(5)=qbpos q(6)=4.d0*qbneg ! CM c---------------------------------------------------------------- c-------------------------------------------------------------------------------------- c*** Distanze interatomiche tra 2 molecole per calcolare parte vdW e elettrostatica index=0 do ii=3,3 do hh=4,5 index=index+1 dist(index)=dsqrt((xxx(ii)-xxx(hh))**2.d0+(yyy(ii)-yyy(hh))**2.d0 < +(zzz(ii)-zzz(hh))**2.d0) enddo enddo CALL EPSIRMATAT(rr2) c epsi1=1.85d0 ! S c rm1=3.770d0 c epsi2=1.08d0 ! C c rm2=3.925d0 c check c write(*,*) epsi1,rm1 c write(*,*) epsi2,rm2 c write(*,*) dipB c stop c Calcola potenziale (vdW) pes1=0.d0 pes2=0.d0 c--------------------------------------------------------------- c coppia S-He epsi=epsi1 rm=rm1 betaS=8.0d0 c if(rr2.ge.rrcs)then c beta=betaS-0.36d0*(rr2-rrcs)**(1.d0/3.d0) c else c beta=betaS-0.36d0*(rrcs-rr2)**(1.d0/3.d0) c endif c new formulation ! beta=betaS+2.d0*(-(rr2-rrcs)+(rr2-rrcs)**2.d0 ! #+2.d0/3.d0*(rr2-rrcs)**3.d0) beta=betaS+2.d0*(-(rr2-rrcs)+0.5d0*(rr2-rrcs)**2.d0 #+0.5d0*(rr2-rrcs)**3.d0 - 0.25d0*(rr2-rrcs)**4.d0) nlj=6 mlj=beta+4.0*(dist(1)/rm)**2 pes1=pes1+epsi*(nlj/(mlj-nlj))*(rm/dist(1))**mlj- #(mlj/(mlj-nlj))*epsi*(rm/dist(1))**nlj c--------------------------------------------------------------- c coppia C-He epsi=epsi2 rm=rm2 betaC=8.2d0 c if(rr2.ge.rrcs)then c beta=betaC-1.55d0*(rr2-rrcs)**(1.d0/3.d0) c else c beta=betaC-1.55d0*(rrcs-rr2)**(1.d0/3.d0) c endif c new formulation ! beta=betaC+6.d0*(-(rr2-rrcs)+(rr2-rrcs)**2.d0 ! #+2.d0/3.d0*(rr2-rrcs)**3.d0) beta=betaC+6.d0*(-(rr2-rrcs)+0.5d0*(rr2-rrcs)**2.d0 #+0.5d0*(rr2-rrcs)**3.d0 - 0.25d0*(rr2-rrcs)**4.d0) nlj=6 mlj=beta+4.0*(dist(2)/rm)**2 pes2=pes2+epsi*(nlj/(mlj-nlj))*(rm/dist(2))**mlj- #(mlj/(mlj-nlj))*epsi*(rm/dist(2))**nlj c--------------------------------------------------------------- c Induzione c Vind =-1389d0*3.d0*(dcos(thetaradb)**2+1.d0)/(R**6) Vind = -(alfaA*DIPB**2)*(1.d0+p2(thetaradb)) #*27211.4d0*0.529177d0**6 #/R**6.d0 c Correzione anisotropia repulsione Vrep=70000d0*dexp(-3.d0*R)*dsin(thetaradb)**6 vtot=pes1+pes2 #+Vind #+Vrep return end c----------------------------------------------------------------------- SUBROUTINE EPSIRMATAT(rrb1) IMPLICIT REAL*8 (A-H,O-Z) COMMON/PARAM2/epsi1,epsi2,rm1,rm2 COMMON/PARAM3/C6medio1,C6medio2,alfamedioA,alfamedioB1,alfamedioB2 #,alfamedioAeff,alfamedioB1eff,alfamedioB2eff c13=1.d0/3.d0 ELMOLA=2.d0 ! per He ELMOLB1=6.8d0 ! per S in CS ELMOLB2=4.2d0 ! per C in CS amediaA=0.2d0 ! polarizzabilità He alfamedioA=amediaA ! per He c----------------------------------------------------------------------- call POLAR(rrB1,aparB,aperB,amediaB,deltaB) !----------------check c write(*,*) amediaB c stop alfamedioB1=amediaB*0.538d0 ! per S alfamedioB2=amediaB*0.462d0 ! per C c---------------------------------------------------------- c primo monomero c He alfamedioAeff=alfamedioA c seconda molecola c CS alfamedioB1eff=alfamedioB1*1.37d0 ! per S in CS c alfamedioB2eff=alfamedioB2*2.303d0 ! per C in CS alfamedioB2eff=alfamedioB2*1.97d0 ! per C in CS c------------------------------------------------------------------ c Rm1 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rm1=(alfamedioAeff**(1.d0/3.d0)+alfamedioB1eff**(1.d0/3.d0))/ #((alfamedioAeff*alfamedioB1eff)**0.095d0) rm1=rm1*1.767d0 c Rm2 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rm2=(alfamedioAeff**(1.d0/3.d0)+alfamedioB2eff**(1.d0/3.d0))/ #((alfamedioAeff*alfamedioB2eff)**0.095d0) rm2=rm2*1.767d0 c correzione per best fit c rm1=rm1/FIT1*1.015d0 c rm2=rm2*FIT2*1.015d0 c---------------------------------------------------------------- c C6 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c6medio1=(alfamedioA*alfamedioB1)/ #(dsqrt(alfamedioA/ELMOLA)+dsqrt(alfamedioB1/ELMOLB1)) c6medio1=c6medio1*15700.d0 c6medio2=(alfamedioA*alfamedioB2)/ #(dsqrt(alfamedioA/ELMOLA)+dsqrt(alfamedioB2/ELMOLB2)) c6medio2=c6medio2*15700.d0 c correzione per best fit c c6medio1=c6medio1/FITC6 c correzione per best fit c c6medio2=c6medio2/FITC6 c epsi ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c 1 c----------------------------------------------------------- epsi1=0.72d0*c6medio1/rm1**6.d0 c 2 c---------------------------------------------------------- epsi2=0.72d0*c6medio2/rm2**6.d0 c correzione per best fit c epsi1=epsi1*FITEPS1/1.10d0 c epsi2=epsi2*FITEPS2/1.10d0 return end c----------------------------------------------------- c========espressioni dei pol. di Legendre=================== REAL*8 FUNCTION P1(x) IMPLICIT REAL*8(A-H,O-Z) p1 = dcos(x) RETURN END REAL*8 FUNCTION P2(x) IMPLICIT REAL*8(A-H,O-Z) p2 = 0.5d0*(3.0d0*dcos(x)*dcos(x) - 1.0d0) RETURN END REAL*8 FUNCTION P3(x) IMPLICIT REAL*8(A-H,O-Z) p3 = 0.5d0*(5.d0*(dcos(x))**3-3.d0*dcos(x)) RETURN END SUBROUTINE POLAR(rr,apar,aper,amedia,delta) c programma calcolo polarizzabilita' CS c in funzione della distanza intermolecolare implicit real*8(a-l,o-z) dimension abond(70,10) c-------------------------------------- c ordine di legame - bo c distanza di legame - bl c alfa atomo 1 - alf1 c alfa atomo 2 - alf2 c elettroni di valenza atomo 1 - ev1 c elettroni di valenza atomo 2 - ev2 c elettroni di non legame atomo 1 - enb1 c elettroni di non legame atomo 2 - enb2 bo=3.0d0 bl=1.5346d0 alf1=1.76d0 alf2=2.90d0 evt1=4.00d0 evt2=6.00d0 enb1=1.00d0 enb2=3.00d0 c-------------------------------------- C COSTANTI c13=0.333333333d0 c12=0.500000000d0 c16=1.0d0/6.0d0 c14=1.0d0/4.0d0 c23=2.0d0/3.0d0 c32=1.d0/c23 C PARAMETRI a1=alf1/evt1 a2=alf2/evt2 at=a1+a2 alfa=alf1**c13+alf2**c13 a13=a1**c13 a23=a2**c13 c------------------------------------- c loop sulla distanza intermolecolare c c rr=0.100d0 c do 111 j=1,100 c------------------------------------- c oppure distanza fissa di equilibrio c rr=bl c------------------------------------- c fattore di contrazione fc=rr/alfa c fattore di spegnimento fs=rr/bl c fattore novita' fn=rr/(a13+a23) c termini esponenziali esp1=2.000d0*fc-1.000d0 c fe1=dexp(-1.000d0*fs*esp1) c 2o cambio fe1=dexp(-1.000d0*(fs-1.d0)*2.d0*fc) esp2=0.7500d0*(fs-1.00d0)*(1.00d0-fc)*(fs**c12)*alfa**c12 fe2=dexp(esp2) esp3=2.00d0*(1.00d0-fs)/bo fe3=dexp(esp3) c esp4=dabs(esp1) c 1o cambio nel valore assoluto esp4 esp4=esp1 fe4=dexp(-0.75d0*esp4*(fs**c12)*alfa**c12) c--------------------------------------- c alfa media legame c--------------------------------------- c fattore campana cex=(fs-1.0d0)*fe2/3.0d0 camp=1.0d0+cex c ordine di legame effettivo c b2=fs+1.00d0 c b2=1.00d0/b2 c b1=dabs(bo-2.0d0) c b1=3.0d0*(bo**b1)/4.0d0 c 3o cambio b2=4.d0*(fc**2)+1.00d0 b2=(fs+2.00d0)/b2 b1=((bo-2.d0)**2)-1.d0/8.d0 b1=1.0d0*(bo**b1)/4.0d0 bof=bo-b1*b2*fe1 c--------------------------------------- c polarizzabilita media (amedia) alfnb=a1*enb1+a2*enb2 alfv=at*bof amedia=(alfv+alfnb)*camp c--------------------------------------- c anisotropia (delta) c--------------------------------------- pre1=c32*alfa**c12*fn**2.d0 pre21=fe3*(enb1+enb2)/2.0d0+(enb1*enb2)**c12 pre21=(pre21/(evt1+evt2)) pre24=bo**((bo-4.d0)/(bo+1.0d0)) pre2=(1.0d0+pre24*pre21)**(3.d0) pre4=2.0d0*fc/(bo**c13)-dabs(alf1-alf2)/(alf1+alf2) delta=fe4*amedia*pre1*pre4/pre2 c--------------------------------------- c parallela--perpendicolare c--------------------------------------- apar=amedia+c23*delta aper=amedia-c13*delta c--------------------------------------- 1 continue return end