1 / 1

jose_morse

marva
Download Presentation

jose_morse

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. From MANAGER@xalfs.mst.lanl.gov Thu Aug 10 12:46:09 1995 To: bridges@cats.UCSC.EDU Subject: code for morse potential subroutine XMORSE (ISHELL,IERR) c purpose to take exact thermal average for single scattering c exafs contribution from a static bond, c using a model potential (morse-potential) c to describe the motion of the bond. c it can be done in the classical limit or using a full c quantum mechanical formalism. c coded by i. batistic and j. mustre (los alamos) c version 1.0 march 1990 c version 2.0 april 1991 include 'EXAFS_FIT.INC' integer ishell,ierr integer kkk parameter (kkk=400) c pi is defined in constants.cmn which is included in exafs_fit.inc c common/const/pi c real pi common/oscil/xx(kkk),dd(kkk),sd(kkk),zz(kkk,kkk),V(kkk) real xx,dd,sd,zz,V common/oscres/mesh,xmm,dx integer mesh real xmm,dx common/param/NH,T integer NH real T common/param1/a,b,amass,cFLAG,alpha,beta real a,b,amass,cFLAG,alpha,beta common/dat/ak(361),hi(361),rLJ(200),sr(200),si(200),sa(200) real ak,hi,rlj,sr,si,sa common/class/xmaxi real xmaxi common/scatt/r0,r00,AMPI,PHSI real r0,r00,ampi,phsi integer mesh1 real fac2,r000,akmax c c pi=3.14159265359 c mesh=0 mesh1=161 NPTS_R = MESH1 c SIGMA(ISHELL) is an extra debye-waller factor to take into account c structural disorder c U0(ISHELL) is the depth of the potential divided by 10 c alpha and beta define the lennard-jones potential c FLAG(ISHELL)=-1 quantum mechanical case, FLAG(ISHELL)=0 classical limit c amass is dimensionless mass related to real mass M by: c M=amass*hbar^2/(e0*r000) c e0 is an overall energy scale choosen such that input: c TEMP(ISHELL)erature t corresponds to real TEMP(ISHELL)erature T/100 in K c e0=kb*T/t c the constant 0.02058 alows the mass used here (amass) to be c entered in atomic units amass=amassP(ISHELL)*0.02058 a=10.*abs(U0(ISHELL)) b=abs(SIGMA(ISHELL)) cFLAG=FLAG(ISHELL) beta=abs(BETAP(ISHELL)) alpha=abs(ALPHAP(ISHELL)) beta=alpha write (*,*) 'morse version 2.0 last modified 8/91 j.m.' write (*,*) 1 'alpha,beta,a,DISTANCE', alpha,beta,U0(ISHELL),DISTANCE(ISHELL) write(*,*)'TEMP(ISHELL),mass=',(tEMP(ISHELL)*100), 1 (amass/0.02058) c fac2 appropriate for morse potential fac2=log(alpha/beta)*(1./(2*alpha-beta)) c r0=abs(DISTANCE(ISHELL)+fac2) c fix length scale r000 to be 0.1 A (arbitrary choice) r000=0.1 r00=r000 c classical limit (NH=0), finite TEMP(ISHELL)erature quantum regime (NH=-1) c QM-case NH state (NH>0) c use parameter c to indicate qm or classical case NH=int(cFLAG) c TEMP(ISHELL)erature entered is assumed to be TEMP(ISHELL) multiplied by 100 t=TEMP(ISHELL) c maximum value of the x-range, fixed later on in oscill1 c and sink subroutines. xmaxi=8.d0 c if (NH .eq. 0) go to 200 102 call schroed_MORSE(a,alpha,beta,amass,mesh1,ISHELL) write (*,*) 'r0, x_max, kmax = ',r0,xmm,pi/xmm 200 continue xmaxi=8.0 201 if (NH.le.0) then 202 if (T.le.0.) T=1. if (NH.eq.0) then write (*,*) 'classical limit' xmaxi=8.0 end if end if c 203 if (r000.eq.0) then r00=0.1/xmm else r00=r000 end if 204 akmax=amin1(4*pi/xmm,pi/dx,6*pi/xmaxi)/(2.*r00) write(*,*) 'NH,TEMP(ISHELL)= ', NH,TEMP(ISHELL)*100. write(*,*) 'a,b,c,xmaxi=',a,b,c,xmaxi write(8,*) 'NH,TEMP(ISHELL)= ', NH,TEMP(ISHELL)*100. 205 call calc(ISHELL) 1000 format(A1) 1100 format(2(A20)) 1200 format(2(E20.5E4)) return End c c23456789012345678901234567890123456789012345678901234567890123456789012 c23456789012345678901234567890123456789012345678901234567890123456789012 c Subroutine schroed_MORSE(a,alpha,beta,amass,mesh1,ISHELL) c purpose to solve one-dimensional Schrodinger equation for c a model potential defined inside include 'EXAFS_FIT.INC' real a,alpha,beta,amass integer mesh1,ishell integer kkk parameter (kkk=361) c common/oscil/xx(kkk),dd(kkk),sd(kkk),zz(kkk,kkk),V(kkk) real xx,dd,sd,zz,v common/oscres/mesh,xmm,dx integer mesh real xmm,dx common/scatt/r0,r00,AMPI,PHSI real r0,r00,ampi,phsi real fac2,xmm1,umin,fac,r1,u integer i,j,k,iflag,ierr c mesh=min0(mesh1,kkk) c fac2=(alpha/(2*beta))**(1./(alpha-beta)) C DISTANCE(ISHELL)=r0*fac2 xmm1=8. xmm=xmm1 Umin=0. c kinetic energy term in Schrodinger equation (hbar^2/2m)* (d psi/dx)^2 c gives rise to a subdiagonal part sd=fac/2 and a diagonal contribution c fac that is added to U, dd=U+fac. dx=2.*xmm1/real(mesh-1.) fac=1./(amass*dx*dx) do 100 k=1,mesh xx(k)=-xmm1+dx*real(k-1) r1=abs(r0+xx(k)*r00) c define potential c morse potential V(k)=a*(exp(-2*alpha*(r1-r0))-2*exp(-(r1-r0)*beta)) Umin=amin1(V(k),Umin) U=V(k) dd(k)=U+fac if (k.lt.mesh) sd(k+1)=-0.5*fac 100 continue if (iflag .ne. 1)then c do 102 k=1,mesh c write(10,101) (0.1*xx(k)),(100*V(k)) c 101 format(1x,2(e12.6,1x)) 102 continue endif sd(1)=0. c do 110 i=1,mesh do 120 j=1,mesh zz(i,j)=0. 120 continue zz(i,i)=1. 110 continue c call IMTQL2(kkk,mesh,dd,sd,zz,ierr) write (*,*) 'IMTQL err = ',ierr write (*,*) 'U0(ISHELL),alpha,beta=',(100*Umin),alpha,beta write (*,*) 'enerqs = ',((100*(dd(i)-Umin)),i=1,15) write (8,*) 'U0(ISHELL),alpha,beta=',(100*Umin),alpha,beta write (8,*) 'enerqs = ',((100*(dd(i)-Umin)),i=1,15) if (WRITE_POT_PARMS) then c U0(ISHELL) = 100. * UMIN ALPHAP(ISHELL) = ALPHA BETAP(ISHELL) = BETA do I=1,15 ENERGY_QS(I,ISHELL) = 100. * (DD(I) - UMIN) enddo endif c iflag=1 return End c23456789012345678901234567890123456789012345678901234567890123456789012 c real Function sin3_MORSE(ak,write_pot_parms,T,Ishell) c classical limit case logical WRITE_POT_PARMS real ak,T integer kkk real theta c parameter (kkk=361, theta=5.58) c c common/scatt/r0,r00,a0,a1,a2,c0,c3 common/scatt/r0,r00,AMPI,PHSI real r0,r00,ampi,phsi common/param1/a,b,amass,cFLAG,alpha,beta real a,b,amass,cFLAG,alpha,beta common/class/xmaxi real xmaxi c fac2=(alpha/(2*beta))**(1./(alpha-beta)) c DISTANCE(ISHELL)=r0*fac2 fi=PHSI sin3_morse=0. psin=0. 101 format(1x,4(e12.6,1x)) 11 do 12 j=1,800 x=(j-361)*xmaxi/real(361.) r1=abs(r0+x*r00) c morse potential v=a*(exp(-2*alpha*(r1-r0))-2*exp(-(r1-r0)*beta)) c boltzmann factor [g(r)] wth=exp(-v/T) c integrate dr chi(r)*g(r) c include in chi a static disorder debye-waller factor sin3_morse=sin3_morse+wth*sin(2.*ak*r1+fi)* 1 exp(-2.*(b*ak)**2)/(r1**2) psin=psin+wth 12 continue c include amplitude that is independent of r sin3_morse=AMPI*sin3_morse/psin iflag=iflag+1 c plot potential +radial distribution function if (iflag .eq. 1 .or. iflag .eq. 303)then sum1=0. sum2=0. sum3=0. sum4=0. pssin=psin 14 do 15 j=1,800 x=(j-361)*xmaxi/real(361.) r1=abs(r0+x*r00) c morse potential v=a*(exp(-2*alpha*(r1-r0))-2*exp(-(r1-r0)*beta)) c boltzmann weight wsth=exp(-v/T) c calculate moments from the rdf(wsth) sum2=sum2+(x**2*wsth/pssin) sum3=sum3+(x**3*wsth/pssin) sum1=sum1+(x*wsth/pssin) sum4=sum4+(x**4*wsth/pssin) if (WRITE_POT_PARMS) then GRID_ADJ = X call WRITE_POTENTIAL3 (J,R00,GRID_ADJ,V,WSTH,PSSIN,ISHELL) C WRITE (5,9000) J,R00,X,V,WSTH,PSSIN C9000 FORMAT (' J R00 X V WSTH PSSIN',I5,5F12.5) endif 15 continue call WRITE_MOMENTS (R00,SUM1,SUM2,SUM3,SUM4,ISHELL) WRITE_POT_PARMS = .false. c 15 write(9,101)(r00*x),(v*100),(500.625*wsth/pssin) c write(9,101)(r00*x),(r00*sum1),(r00*sqrt(sum2)), c 1 (r00*(abs(sum3))**1/3.),(r00*sum4**0.25) c write(*,*)' delta r, delta r2, delta r3, delta r4',(r00*sum1), c 1 (r00*sqrt(sum2)), c 2 (r00*(abs(sum3))**1/3.),(r00*sum4**0.25) c write(8,*)' delta r, delta r2, delta r3, delta r4',(r00*sum1), c 1 (r00*sqrt(sum2)), c 2 (r00*(abs(sum3))**1/3.),(r00*sum4**0.25) else return endif return End

More Related