CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PROGRAM E1E12 C Bohr-Sommerfeld quantization for bound states of the C Lennard-Jones, or Morse, Potential CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL ARCHON !search for bound states END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE ARCHON C finds the bound states of the Lennard-Jones, or Morse, potential C from the Bohr-Sommerfeld quantization rule CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Local variables: REAL E1 !current value of energy REAL X1,X2 !current turning points REAL F1 !f=action/2 - pi/2 -ilevel*pi INTEGER ILEVEL !current level C Global variables: INCLUDE 'PAR2.E1' C REAL ENERGY(0:MAXLVL) !energy of bound state REAL XIN(0:MAXLVL) !inner turning point REAL XOUT(0:MAXLVL) !outer turning point C Global variables: INCLUDE 'PAR12.E1' PRINT*,' IOPT = 0 (1) for Lennard-Jones (Morse) potential' PRINT*,' enter value of IOPT' READ*,IOPT IF(IOPT.EQ.0) THEN POTMIN=1.12246 ELSE POTMIN=0.74166 ENDIF IF(IOPT.EQ.1) THEN PRINT*,' enter value of BETA' READ*,BETA ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C OPEN(11, FILE='LIXO.OUT',STATUS='UNKNOWN') C search for bound states E1=-1. !begin at the well bottom F1=-PI/2 !the action is zero there C C find the NLEVEL bound states DO 100 ILEVEL=0,NLEVEL-1 C CALL SEARCH(ILEVEL,E1,F1,X1,X2) !search for eigenvalue ENERGY(ILEVEL)=E1 !store values for this state XIN(ILEVEL)=X1 XOUT(ILEVEL)=X2 C C C output WRITE(11,*)'N= ',ILEVEL,' ENERGY= ',VH2*ENERGY(ILEVEL) PRINT*,'N= ',ILEVEL,' ENERGY= ',VH2*ENERGY(ILEVEL) F1=F1-PI !guess to begin search for next level C 100 CONTINUE C C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SEARCH(N,E1,F1,X1,X2) C finds the N'th bound state C E1 is passed in as initial guess for the bound state energy C and returned as the true bound state energy with turning points C X1 and X2 C F1 is the function which goes to zero at a bound state C F1 = action/2-(n+1/2)*pi CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Input/Output variables: INTEGER N !current level (input) REAL E1,E2 !trial energies (I/O) REAL F1,F2 !f=action/2-pi*(n+1/2) (I/O) REAL S !action (output) REAL X1,X2 !turning points (output) C Local variables: REAL DE !increment in energy search C Global variables: INCLUDE 'PAR12.E1' CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C guess the next energy in order to begin search E2=E1+ABS(E1)/4. DE=2*ETOL C C use secant search to find the bound state 50 IF (ABS(DE) .GT. ETOL) THEN CALL ACTION(E2,X1,X2,S) !S at new energy F2=S-(N+.5)*PI !F at new energy IF (F1 .NE. F2) THEN !calculate new DE DE=-F2*(E2-E1)/(F2-F1) ELSE DE=0. END IF C E1=E2 !roll values F1=F2 E2=E1+DE !increment energy IF (E2 .GE. 0) E2=-ETOL !keep energy negative GOTO 50 END IF C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE ACTION(E,X1,X2,S) C calculates the (action integral)/2 (S) and the classical turning C points (X1,X2) for a given energy (E) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Input/Output variables: REAL E !energy (input) REAL S !action (output) REAL X1,X2 !turning points (output) C Local variables: REAL DX !increment in turning point search REAL H !quadrature step size REAL SUM !sum for integral INTEGER IFAC !coefficient for Simpson's rule INTEGER IX !index on X REAL X !current X value in sum REAL POT !potential as a function of X (function) C Global variables: INCLUDE 'PAR12.E1' CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C find inner turning point; begin search at the well bottom X1=POTMIN DX=.1 50 IF (DX .GT. XTOL) THEN X1=X1-DX !use simple search, going inward IF (POT(X1) .GE. E) THEN X1=X1+DX DX=DX/2 END IF GOTO 50 END IF C C find the outer turning point; begin search at the well bottom X2=POTMIN DX=.1 120 IF (DX .GT. XTOL) THEN X2=X2+DX !use simple search going outward IF (POT(X2) .GE. E) THEN X2=X2-DX DX=DX/2 END IF GOTO 120 END IF C C Simpson's rule from X1+H to X2-H IF (MOD(NPTS,2) .EQ. 1) NPTS=NPTS+1 !NPTS must be even H=(X2-X1)/NPTS !step size SUM=SQRT(E-POT(X1+H)) IFAC=2 DO 200 IX=2,NPTS-2 X=X1+IX*H IFAC=6-IFAC !alternate factors SUM=SUM+IFAC*SQRT(E-POT(X)) 200 CONTINUE SUM=SUM+SQRT(E-POT(X2-H)) SUM=SUM*H/3 C C special handling for sqrt behavior of first and last intervals SUM=SUM+SQRT(E-POT(X1+H))*2*H/3 SUM=SUM+SQRT(E-POT(X2-H))*2*H/3 S=SUM*GAMMA C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL FUNCTION POT(X) C evaluates the Lennard-Jones potential at X CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Passed variables: REAL X C Global variables: INCLUDE 'PAR12.E1' CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C If you change the potential, normalize to a minimum of -1 C and change the value of POTMIN in input file PAR12.E1 to the C new equilibrium position (i.e. the X value at which the force is zero) C Lennard-Jones potential in scaled variables IF(IOPT.EQ.0) THEN POT=4*(X**(-12)-X**(-6)) ELSE POT=(1.-EXP(-BETA*(X-POTMIN)))**2-1. ENDIF RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC