program anharm
c   Semiclassical expansion for arbitrary potential v0 + v2 x^2 + v3 x^3 + ...
      implicit real*8(a-h,o-z)
      dimension f(999),e(999),x(3333),t(710,710)
c   N <~200, n <~100
      one=1
      two=2
      open(11,FILE='../../system/temp/anharm.inp')
      read(11,*)nqx
      read(11,*)morder
c Zero order - minimum of the potential
      read(11,*)v0
      e0=v0
c The first derivative, should be zero
      read(11,*)v1
c 1/2 of the second derivative, should be positive or zero (see below)
      read(11,*)v2
c If v2=0, then set omega=1 and print only RSPT series (starting with (nqx+1/2)omega)
      omega=1
      if(v2.ne.0)omega=Sqrt(2*v2)
c Reading anharmonicity and scaling the potential to 1/2 x^2 + f_1 x^3 + f_2 x^4 + ...
      nm=morder-1
      nmd=2*nm
      do 21 i=1,nmd
      read(11,*)vi
  21  f(i)=vi/omega**((i+2)/two)
      close(11)
c First order - quantum harmonic vibrations      
      e1=(nqx+one/2)*omega
c Printing E_0, E_1
      open(12,FILE='../../system/temp/anharm.out')
c For quadruple precision, change everywhere es23.15 to es39.31
 100  format('E<SUB>0</SUB>   =',es23.15,' (Minimum of the potential)')
      if(v2.ne.0)print100,e0
      write(12,*)e0
 110  format('E<SUB>1</SUB>   =',es23.15,
     - ' (Quantum harmonic-oscillator energy)')
      if(v2.ne.0)print110,e1
 210  format('E<SUB>0</SUB>   =',es23.15,
     - '  (Unperturbed harmonic-oscillator energy)')
      if(v2.eq.0)print210,e1
      write(12,*)e1
      na=nqx+1
      SQ2=1/sqrt(TWO)
      NP=MIN0(3*NM+1,NA)
      MM=NP+3*NM
      MM2=MM+2
      MMH=(MM+1)/2
      NPT=(NP+1)/2
      NF=1
      DO4N=1,NMD
      DO18I=1,MMH
   18 T(N,I)=0
      N9=N-1
      IF(F(N).NE.0.d0)NF=N
      IN=1
      IF((NP+N)/2*2.EQ.NP+N)IN=2
      MA1=MIN0(NP+3*N,NP+3*(NMD-N))
      MI1=MAX0(NP-3*N,NP-3*(NMD-N),IN)
      DO5NQ=MI1,MA1,2
      A=0
      DO17I=1,MM2
  17  X(I)=0
      NQT=(NQ+1)/2
      NB=NA+NQ-NP
      IF(NQ.GT.1)X(NQ-1)=(NB-1)*(NB-2)
      X(NQ+1)=2*NB-1
      X(NQ+3)=1
      IF(N.EQ.1)GOTO12
      DO6NX=1,N9
      NT=N-NX
      IF(NX.GT.NF)GOTO11
      IN=1
      IF((NP+NT)/2*2.EQ.NP+NT)IN=2
      MA2=MIN0(NP+3*NT,NQ+NX+2)
      MI2=MAX0(NP-3*NT,NQ-NX-2,IN)
      IF(MI2.GT.MA2)GOTO12
      C=0
      DO9NR=MI2,MA2,2
      NRT=(NR+1)/2
      B=X(NR)+(NR-NP+NA)*X(NR+2)
      X(NR+1)=B
  9   C=C+B*T(NT,NRT)
      A=A-C*F(NX)/2*SQ2**NX
  11  NXH=NX/2
      IF(NXH*2.NE.NX)GOTO6
      A=A+E(NXH)*T(NT,NQT)
   6  CONTINUE
  12  A=A-F(N)/2*SQ2**N*(X(NP)+NA*X(NP+2))
      IF(NQ.NE.NP)A=A/(NQ-NP)
      T(N,NQT)=A
   5  CONTINUE
      NH=N/2
      IF(NH*2.NE.N)GOTO4
      E(NH)=-T(N,NPT)
  10  format('E<SUB>',i1,'</SUB>   =',es23.15)
  20  format('E<SUB>',i2,'</SUB>  =',es23.15)
  30  format('E<SUB>',i3,'</SUB> =',es23.15)
      nh1=nh
      if(v2.ne.0)nh1=nh+1
      if(nh1.le.9)print10,nh1,omega*e(nh)
      if(nh1.ge.10 .and. nh1.le.99)print20,nh1,omega*e(nh)
      if(nh1.ge.100)print30,nh1,omega*e(nh)
      write(12,*)omega*e(nh)
      T(N,NPT)=0
   4  CONTINUE
      close(12)
      stop
      END