公開プログラム(第0.3版、3/11、1998)

やさしいバンド計算プログラムの作り方へ[戻る][旧版][ガイドへ][Top][Forbidden redistribute](再配布禁止)[無保証]

C     @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C     @@ BAND STRUCTURE AND TOTAL ENERGY CALCULATION @@
C     @@           (ORIGINAL FROM  DR. H. ISHIDA)    @@
C     @@ 1989  5/4 (START) BY KOBAYASHI KAZUAKI      @@
C     @@            MODIFIED BY Y. MORIKAWA(SYMMETRY)@@
C     @@ 1992  7/14 MODIFIED TO STRESS CALCULATION   @@
C     @@ 1993  3/16 MODIFIED TO CONSTANT PRESSURE    @@
C     @@            AUTOMATIC OPT (ANDERSEN ARG)     @@
C     @@ 1994 10/25 MODIFIED TO INTRODUCTION OF PCC  @@
C     @@ 1994 11/22 MODIFIED TO CONSTANT PRESSURE    @@
C     @@            AUTOMATIC OPT (PARRINELLO-RAHMAN)@@
C     @@ 1994 12/16 MODIFIED TO INPUT/OUTPUT OF RECAL@@
C     @@ 1995  1/2  MODIFIED TO VPP-PARALEEL VERSION @@
C     @@ 1995  2/13             VPP-PARALLEL VER 1.01@@
C     @@ 1995  6/2              VPP-PARALLEL VER 1.1 @@
C     @@ 1995  7/17 MODIFIED TO ADD NON-LOCAL D  1.2 @@
C     @@ 1997  5/1  NEW VERSION VER 1.2A.01 (START)  @@
C     @@ 1997  5/21 MAIN03 ---> MAIN04 (1.2A.04)     @@
C     @@ 1998  3/11 ADD SUB DNV,BUBBLE       VER 1.2B@@
C     @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
      IMPLICIT REAL(A-H,O-Y)
      IMPLICIT COMPLEX(Z)
      INCLUDE 'PACVPP'
C*****FOR PCC
      INTEGER   KTPCC(KTYP)
      DIMENSION RHPC(KIZA,KTYP)
C*****FOR KPMIRB
      DIMENSION QWGT(KNV3)
C*****FOR SYMM
      INTEGER   IRX(125),IRY(125),IRZ(125)
      DIMENSION RX(125),RY(125),RZ(125),RR(125)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      INCLUDE 'REINIT'
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
!XOCL PARALLEL REGION
      KOPR=KO
      WRITE (6,*) 'START OF PROGRAM (REVPE_D.F)'
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL SIMP
      CALL TIME(ITIME1,'SIMP  ',2)
C^^^^ READ UNIT = 1 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL INPUT1(ALFA)
      CALL TIME(ITIME1,'INPUT1',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL INFOUT(ICAR,ICAF,IPRE,IPCC,IMD,IOVE,IRES,IORW,ISSS
     &     ,ISTMD,ICOST,IOUT,IMDI,NBAND1,NBAND2
     &     ,IREC8,IREC9,KREC8,KREC9,KBZTYP
     &     ,WIDTH,CMIX,DTIMEL,DTIMIO,DTIMUC,KTPCC)
      CALL TIME(ITIME1,'INFOUT',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      ITER=0
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL LATTIC(RX,RY,RZ,RR,IRX,IRY,IRZ)
      CALL TIME(ITIME1,'LATTIC',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL GSTEP1(IPRE)
      CALL TIME(ITIME1,'GSTEP1',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL BUBBLE(KG,IG1,IG2,IG3,GX,GY,GZ,GR)
C      CALL HPSORT(KG,IG1,IG2,IG3,GX,GY,GZ,GR)
      CALL TIME(ITIME1,'BUBBLG',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      IF(KBZTYP.GE.2) THEN
        CALL TIME(ITIME1,'      ',1)
        CALL SYMM(KOPR,NOPR,KBZTYP,125,6,RX,RY,RZ)
        CALL TIME(ITIME1,'SYMM  ',2)
      END IF
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL GSTEPF
      CALL TIME(ITIME1,'GSTEPF',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL KSTEP
     >     (KNV3,KBZTYP,ALTV,RLTV,
     >      KNVX,KNVY,KNVZ,KNVX2,KNVY2,KNVZ2,6,
     <      KV3,VX,VY,VZ,QWGT)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL BASNUM(QWGT)
      CALL TIME(ITIME1,'BASNUM',2)
C^^^^ PRE-CALCULATION IPRE=1 ^^^^^^^^^^^^^^^^^^^^^
      IF (IPRE.EQ.1) THEN
          WRITE (6,*) '[THIS IS PRE-CALCULATION PROGRAM STOPS]'
          STOP
      END IF
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL FORM
      CALL TIME(ITIME1,'FORM  ',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      IF (IPCC.EQ.1) THEN
          CALL TIME(ITIME1,'      ',1)
          CALL PCC(KTPCC,RHPC,SCHGPC,EPC,PCM)
          CALL TIME(ITIME1,'PCC   ',2)
      ELSE IF (ITER.EQ.0) THEN
          EPC     = 0.D0
          DO 1010 I=1,KNG
            ZRHPC(I)=DCMPLX(0.D0,0.D0)
            ZRRPC(I)=DCMPLX(0.D0,0.D0)
 1010     CONTINUE
      END IF
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL PSEUDO(TOTCH,ETOT1)
CCC      CALL PSELMD
      CALL TIME(ITIME1,'PSEUDO',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL EWVEC(ETOT2)
CCC      CALL EWVMD(ETOT2)
      CALL TIME(ITIME1,'EWALDV',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL INTCHG(ALFA)
      CALL TIME(ITIME1,'INTCHG',2)
C^^^^ LOOP ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 1000 CONTINUE
      ITER=ITER+1
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL XCFFT(1,SCHGPC)
      CALL TIME(ITIME1,'XCFFT ',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^DEBUG^^^^^^^^^^^^^^^^^^^
      IF (ITER.EQ.1) THEN
          CALL TIME(ITIME1,'      ',1)
          CALL KBMAT
          CALL DIAGON(IREC8)
          CALL KBINT
          CALL TIME(ITIME1,'ITER1 ',2)
      ELSE
          CALL TIME(ITIME1,'      ',1)
          CALL MSD(IREC8,IREC9,IMD,ICAR,ISTR)
          CALL TIME(ITIME1,'MSD   ',2)
      END IF
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL FERMI(TOTCH,QWGT,WIDTH)
      CALL TIME(ITIME1,'FERMI ',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL FORZFB
      CALL TIME(ITIME1,'FORCE ',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL CHAVER(IREC8,IREC9,KOPR,NOPR,KBZTYP)
      CALL TIME(ITIME1,'CHAVER',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      IF (IPCC.EQ.1) THEN
          CALL TIME(ITIME1,'      ',1)
          CALL XCFFT(3,SCHGPC)
          CALL TIME(ITIME1,'XCFFT ',2)
      END IF
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL ENERGY(TOTCH,ETOT1,ETOT2,ETONEW,SCHGPC,EPC)
      CALL TIME(ITIME1,'ENERGY',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      CALL TIME(ITIME1,'      ',1)
      CALL CONV2(ETONEW,ETOOLD,IMD,ITERMD,CMIX)
      CALL TIME(ITIME1,'CONV2 ',2)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      IF (ITER.GT.IMD) GO TO 9999
      GO TO 1000
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
!XOCL END PARALLEL
 9999 CONTINUE
      STOP
      END
C---------------------------------------------------------              
      SUBROUTINE TIME(ITIME1,NAME,ISE)                                  
      IMPLICIT REAL(A-H,O-Y)
      IMPLICIT COMPLEX(Z)
      COMMON/IT/     ITER                                               
      CHARACTER NAME*6                                                  
      IF (ITER.LE.2 .OR. ITER.EQ.10) THEN                              
        IF (ISE.EQ.2) THEN                                              
C            CALL CLOCKM(ITIME)                                          
            T=(ITIME-ITIME1)/1000.0D0                                   
            WRITE (6,600) NAME,T                                        
          ELSE                                                          
C            CALL CLOCKM(ITIME1)                                         
             T = 0.0D0
        END IF                                                          
      END IF                                                            
      RETURN                                                            
  600 FORMAT(1H ,'SUB-NAME = ',A10,'  ----------------- TIME = ',F12.6) 
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                      
      SUBROUTINE INFOUT(ICAR,ICAF,IPRE,IPCC,IMD,IOVE,IRES,IORW,ISSS
     &     ,ISTMD,ICOST,IOUT,IMDI,NBAND1,NBAND2
     &     ,IREC8,IREC9,KREC8,KREC9,KBZTYP
     &     ,WIDTH,CMIX,DTIMEL,DTIMIO,DTIMUC,KTPCC)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                      
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      INTEGER   KTPCC(KTYP),KMDTP(KATM)
C                                                                       
      WRITE (6,*) 'IPRE      = ',IPRE
      WRITE (6,*) 'ICAR      = ',ICAR
      WRITE (6,*) 'ICAF      = ',ICAF
      WRITE (6,*) 'IPRI      = ',IPRI
      WRITE (6,*) 'ICONT     = ',ICONT
      WRITE (6,*) 'IPCC      = ',IPCC
      WRITE (6,*) 'KBZTYP    = ',KBZTYP
      WRITE (6,*) 'IMD       = ',IMD
      WRITE (6,*) 'IOVE      = ',IOVE
      WRITE (6,*) 'IMDI      = ',IMDI
      WRITE (6,*) 'IRES      = ',IRES
      WRITE (6,*) 'IORW      = ',IORW
      WRITE (6,*) 'ISSS      = ',ISSS
      WRITE (6,*) 'ISTOP     = ',ISTOP
      WRITE (6,*) 'ISTMD     = ',ISTMD
      WRITE (6,*) 'ICOST     = ',ICOST
      WRITE (6,*) 'IOUT      = ',IOUT
      WRITE (6,*) 'IREC8     = ',IREC8                                  
      WRITE (6,*) 'IREC9     = ',IREC9                                  
      WRITE (6,*) 'KREC8     = ',KREC8                                  
      WRITE (6,*) 'KREC9     = ',KREC9                                  
      WRITE (6,*) 'UP-BAND   = ',NBAND1                                 
      WRITE (6,*) 'DOWN-BAND = ',NBAND2                                 
C     OUTPUT INCLUDE FILE INFORMATION                                   
      WRITE (6,*) 'OUTPUT INCLUDE FILE INFORMATION'                     
      WRITE (6,*) 'KNX,Y,Z   =',KNX,KNY,KNZ                             
      WRITE (6,*) 'KNG,1,11  =',KNG,KNG1,KNG11                          
      WRITE (6,*) 'KNBMX     =',KNBMX
      WRITE (6,*) 'KNVX,Y,Z  =',KNVX,KNVY,KNVZ                          
      WRITE (6,*) 'KVX2,Y2,Z2=',KNVX2,KNVY2,KNVZ2
      WRITE (6,*) 'KEG,KBD1,2=',KEG,KBD1,KBD2                           
      WRITE (6,*) 'KTYP,KATM =',KTYP,KATM                               
      WRITE (6,*) 'EWALD     =',KLX,KLY,KLZ,KMX,KMY,KMZ
      WRITE (6,*) 'DTIMEL    =',DTIMEL
      WRITE (6,*) 'DTIMIO    =',DTIMIO
      WRITE (6,*) 'DTIMUC    =',DTIMUC
      WRITE (6,*) 'WIDTH     =',WIDTH
      WRITE (6,*) 'CMIX      =',CMIX
      WRITE (6,*) 'TEMPSB    =',TEMPSB
      DO 1000 ITYP=1,KTYP
      WRITE (6,*) 'KTPCC     =',KTPCC(ITYP)
 1000 CONTINUE
      DO 1010 I=1,6
      WRITE (6,*) 'PEX = ',I,PEX(I)
 1010 CONTINUE
      DO 2000 I=1,KATM                                                  
      WRITE (6,*) 'AMION = ',AMION(I)
 2000 CONTINUE                                                          
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                      
      SUBROUTINE SIMP                                                   
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                      
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      PARAMETER(KIZA=421)                                               
      COMMON/OM/ OMO(KIZA)                                              
      OMO(1)     = 1.0D0                                                
      DO 1101 N  = 2,KIZA-1,2                                           
      OMO(N)     = 4.0D0                                                
      OMO(N+1)   = 2.0D0                                                
 1101 CONTINUE                                                          
      OMO(KIZA)  = 1.0D0                                                
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                   
      SUBROUTINE INPUT1(ALFA)                                           
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                   
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
C'''''READING THE DATA (UNIT=1)''''''''''''''''''''''''''''''''''       
      REWIND 1                                                          
      READ(1,*) ITEMAX,PPMIX,CONV,GMAX,ICONT                            
      WRITE(6,99) ITEMAX,PPMIX,CONV,GMAX,ICONT                          
   99 FORMAT(' ',I6,F8.4,D12.4,F8.4,I4)                                 
      DO 100 I=1,3                                                      
      READ(1,*) ALTV(1,I),ALTV(2,I),ALTV(3,I)                           
      WRITE(6,200) ALTV(1,I),ALTV(2,I),ALTV(3,I)                        
  100 CONTINUE                                                          
  200 FORMAT(3(F20.10))                                                 
      READ(1,*) KCOTYP                                                  
      WRITE(6,*)KCOTYP,' COORDINATES  0:NORMALIZED  1:CARTESIAN '       
      IF (KCOTYP.EQ.0) THEN                                             
         DO 300 IA=1,KATM                                               
           READ(1,*)    PPOS(IA),QPOS(IA),RPOS(IA)                      
           WRITE(6,200) PPOS(IA),QPOS(IA),RPOS(IA)                      
  300    CONTINUE                                                       
         DO 400 IA=1,KATM                                               
           CATX(IA)=PPOS(IA)*ALTV(1,1)+QPOS(IA)*ALTV(1,2)               
     &                                +RPOS(IA)*ALTV(1,3)               
           CATY(IA)=PPOS(IA)*ALTV(2,1)+QPOS(IA)*ALTV(2,2)               
     &                                +RPOS(IA)*ALTV(2,3)               
           CATZ(IA)=PPOS(IA)*ALTV(3,1)+QPOS(IA)*ALTV(3,2)               
     &                                +RPOS(IA)*ALTV(3,3)               
  400    CONTINUE                                                       
       ELSE IF(KCOTYP.EQ.1) THEN                                        
         DO 500 IA=1,KATM                                               
           READ(1,*)    CATX(IA),CATY(IA),CATZ(IA)                      
           WRITE(6,200) CATX(IA),CATY(IA),CATZ(IA)                      
  500    CONTINUE                                                       
       ELSE                                                             
         WRITE(6,*) ' KCOTYP = ',KCOTYP                                 
         STOP                                                           
      END IF                                                            
      NNATM=0                                                           
      DO 600 IT=1,KTYP                                                  
        READ(1,*)   IATOM(IT),NLSPD(IT)
        WRITE(6,*) 'IATOM,NLSPD=',IATOM(IT),NLSPD(IT)
        READ(1,*)       AICH,ALFA                                       
        WRITE(6,700) IT,AICH,ALFA                                       
        READ(1,*) ACHG(IT),AC(IT,1),AC(IT,2),BC(IT,1),BC(IT,2)          
        WRITE(6,701) ACHG(IT),AC(IT,1),AC(IT,2),BC(IT,1),BC(IT,2)          
        DO 800 IA=NNATM+1,NNATM+IATOM(IT)                               
          AICHG(IA)   = AICH                                            
          KFTYPE(IA)  = IT                                              
  800   CONTINUE                                                        
        NNATM=NNATM+IATOM(IT)                                           
  600 CONTINUE                                                          
  700 FORMAT(' ','TYPE=',I3,' CHARGE=',F8.4,'GAUSS COEFF.=',F8.4)       
  701 FORMAT(' ','AC,BC=',5F12.6) 
      IF (NNATM.NE.KATM) THEN                                           
          WRITE(6,*) ' NNATM= ',NNATM                                   
          STOP                                                          
      END IF                                                            
      KV3  =  KNV3                                                      
C'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''        
      PINIT=PPMIX                                                       
      DO 1 I=1,ITEMAX                                                   
        PMIX(I)=PPMIX                                                   
    1 CONTINUE                                                          
      DO 2 K=1,KV3                                                      
        DO 3 I=1,KEG                                                    
          OCCUP(I,K)=1.0D0                                              
    3   CONTINUE                                                        
    2 CONTINUE                                                          
      DO 4 I=1,KATM                                                     
        VEL(I,1) =0.0D0                                                 
        VEL(I,2) =0.0D0                                                 
        VEL(I,3) =0.0D0                                                 
 4    CONTINUE                                                          
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                
      SUBROUTINE LATTIC                                                 
     <     (RX,RY,RZ,RR,IRX,IRY,IRZ)                                    
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      REAL   RX(125),RY(125),RZ(125),RR(125)                          
      INTEGER  IRX(125),IRY(125),IRZ(125)                               
      NEIBR  = 2                                                        
      NEIBRD = 125                                                      
C                                                                       
      ALEN1=SQRT(ALTV(1,1)**2+ALTV(2,1)**2+ALTV(3,1)**2)                
      ALEN2=SQRT(ALTV(1,2)**2+ALTV(2,2)**2+ALTV(3,2)**2)                
      ALEN3=SQRT(ALTV(1,3)**2+ALTV(2,3)**2+ALTV(3,3)**2)                
      ALMAX=MAX(ALEN1,ALEN2,ALEN3)                                      
      ALF=1.0D0/ALMAX                                                   
      WRITE(6,*) '      <<>>'                
       FFF=    ALTV(1,1)*(ALTV(2,2)*ALTV(3,3)-ALTV(3,2)*ALTV(2,3))      
     &      +  ALTV(2,1)*(ALTV(3,2)*ALTV(1,3)-ALTV(1,2)*ALTV(3,3))      
     &      +  ALTV(3,1)*(ALTV(1,2)*ALTV(2,3)-ALTV(2,2)*ALTV(1,3))      
                          UNIVOL=ABS( FFF )                             
                          WRITE(6,708) UNIVOL                           
708                       FORMAT(1H ,'VOLUME OF A UNIT CELL=',F15.6)    
        FFF=2.D0*PAI/FFF                                                
      RLTV(1,1)=(ALTV(2,2)*ALTV(3,3)-ALTV(3,2)*ALTV(2,3))*FFF           
      RLTV(2,1)=(ALTV(3,2)*ALTV(1,3)-ALTV(1,2)*ALTV(3,3))*FFF           
      RLTV(3,1)=(ALTV(1,2)*ALTV(2,3)-ALTV(2,2)*ALTV(1,3))*FFF           
          RLTV(1,2)=(ALTV(2,3)*ALTV(3,1)-ALTV(3,3)*ALTV(2,1))*FFF       
          RLTV(2,2)=(ALTV(3,3)*ALTV(1,1)-ALTV(1,3)*ALTV(3,1))*FFF       
          RLTV(3,2)=(ALTV(1,3)*ALTV(2,1)-ALTV(2,3)*ALTV(1,1))*FFF       
              RLTV(1,3)=(ALTV(2,1)*ALTV(3,2)-ALTV(3,1)*ALTV(2,2))*FFF   
              RLTV(2,3)=(ALTV(3,1)*ALTV(1,2)-ALTV(1,1)*ALTV(3,2))*FFF   
              RLTV(3,3)=(ALTV(1,1)*ALTV(2,2)-ALTV(2,1)*ALTV(1,2))*FFF   
              WRITE(6,707) ((RLTV(I,J),I=1,3),J=1,3)                    
  707 FORMAT((1H ,3(F10.6,5X)))                                         
              DO 800 I=1,3                                              
              DO 810 J=1,3                                              
                ALINV(I,J)=RLTV(J,I)/(2.D0*PAI)                         
 810          CONTINUE                                                  
 800          CONTINUE                                                  
              DO 703 I=1,3                                              
              DO 703 J=1,3                                              
        FFF = ALTV(I,1)*ALINV(1,J) + ALTV(I,2)*ALINV(2,J)               
     &                             + ALTV(I,3)*ALINV(3,J)               
         WRITE(6,702) I,J,FFF                                           
 702     FORMAT(1H ,'ALTV(',I3,')*ALINV(',I3,')=',F20.15)               
 703     CONTINUE                                                       
C--*--GENERATE TRANSLATIONAL VECTOR                                     
      MM=0                                                              
      DO 500 I=-NEIBR,NEIBR                                             
        DO 510 J=-NEIBR,NEIBR                                           
          DO 520 K=-NEIBR,NEIBR                                         
            MM = MM + 1                                                 
            IRX(MM) = I                                                 
            IRY(MM) = I                                                 
            IRZ(MM) = I                                                 
                                      F1  = DFLOAT(I)                   
                                      F2  = DFLOAT(J)                   
                                      F3  = DFLOAT(K)                   
            RX(MM)  = ALTV(1,1)*F1 + ALTV(1,2)*F2 + ALTV(1,3)*F3        
            RY(MM)  = ALTV(2,1)*F1 + ALTV(2,2)*F2 + ALTV(2,3)*F3        
            RZ(MM)  = ALTV(3,1)*F1 + ALTV(3,2)*F2 + ALTV(3,3)*F3        
            RR(MM)  = SQRT( RX(MM)**2 + RY(MM)**2 + RZ(MM)**2 )         
  520     CONTINUE                                                      
  510   CONTINUE                                                        
  500 CONTINUE                                                          
      IF(MM.NE.NEIBRD) THEN                                             
        WRITE(6,*) ' MM.NE.NEIBRD MM,NEIBRD =',MM,NEIBRD                
      END IF                                                            
      CALL BUBBLE(NEIBRD,IRX,IRY,IRZ,RX,RY,RZ,RR)                       
C      CALL HPSORT(NEIBRD,IRX,IRY,IRZ,RX,RY,RZ,RR)                       
      IF(IPRI.GE.2) THEN                                                
          WRITE(6,*) ' TRANSLATIONAL VECTOR'                            
          WRITE(6,600) (I,IRX(I),IRY(I),IRZ(I),                         
     &                        RX(I) ,RY(I) ,RZ(I) ,RR(I),I=1,NEIBRD)    
      END IF                                                            
  600 FORMAT((' ',4I5,4F12.6))                                          
      WRITE (6,*) 'UNIVOL(IN LATTIC) = ',UNIVOL                         
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                
      SUBROUTINE LATSCA                                                 
     <     (RX,RY,RZ,RR,IRX,IRY,IRZ)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      REAL   RX(125),RY(125),RZ(125),RR(125)                          
      INTEGER  IRX(125),IRY(125),IRZ(125)                               
      NEIBR  = 2                                                        
      NEIBRD = 125                                                      
C                                                                       
      ALEN1=SQRT(ALTV(1,1)**2+ALTV(2,1)**2+ALTV(3,1)**2)                
      ALEN2=SQRT(ALTV(1,2)**2+ALTV(2,2)**2+ALTV(3,2)**2)                
      ALEN3=SQRT(ALTV(1,3)**2+ALTV(2,3)**2+ALTV(3,3)**2)                
      ALMAX=MAX(ALEN1,ALEN2,ALEN3)                                      
      ALF=1.0D0/ALMAX                                                   
C     WRITE(6,*) '      <<>>'                
       FFF=    ALTV(1,1)*(ALTV(2,2)*ALTV(3,3)-ALTV(3,2)*ALTV(2,3))      
     &      +  ALTV(2,1)*(ALTV(3,2)*ALTV(1,3)-ALTV(1,2)*ALTV(3,3))      
     &      +  ALTV(3,1)*(ALTV(1,2)*ALTV(2,3)-ALTV(2,2)*ALTV(1,3))      
                          UNIVOL=ABS( FFF )                             
                          WRITE(6,708) UNIVOL                           
  708                     FORMAT(1H ,'VOLUME AT LATSCA CELL=',F15.6)    
        FFF=2.D0*PAI/FFF                                                
      RLTV(1,1)=(ALTV(2,2)*ALTV(3,3)-ALTV(3,2)*ALTV(2,3))*FFF           
      RLTV(2,1)=(ALTV(3,2)*ALTV(1,3)-ALTV(1,2)*ALTV(3,3))*FFF           
      RLTV(3,1)=(ALTV(1,2)*ALTV(2,3)-ALTV(2,2)*ALTV(1,3))*FFF           
          RLTV(1,2)=(ALTV(2,3)*ALTV(3,1)-ALTV(3,3)*ALTV(2,1))*FFF       
          RLTV(2,2)=(ALTV(3,3)*ALTV(1,1)-ALTV(1,3)*ALTV(3,1))*FFF       
          RLTV(3,2)=(ALTV(1,3)*ALTV(2,1)-ALTV(2,3)*ALTV(1,1))*FFF       
              RLTV(1,3)=(ALTV(2,1)*ALTV(3,2)-ALTV(3,1)*ALTV(2,2))*FFF   
              RLTV(2,3)=(ALTV(3,1)*ALTV(1,2)-ALTV(1,1)*ALTV(3,2))*FFF   
              RLTV(3,3)=(ALTV(1,1)*ALTV(2,2)-ALTV(2,1)*ALTV(1,2))*FFF   
              WRITE(6,707) ((RLTV(I,J),I=1,3),J=1,3)                    
  707 FORMAT((1H ,3(F10.6,5X)))                                         
              DO 800 I=1,3                                              
              DO 810 J=1,3                                              
                ALINV(I,J)=RLTV(J,I)/(2.D0*PAI)                         
  810         CONTINUE                                                  
  800         CONTINUE                                                  
              DO 703 I=1,3                                              
              DO 703 J=1,3                                              
        FFF = ALTV(I,1)*ALINV(1,J) + ALTV(I,2)*ALINV(2,J)               
     &                             + ALTV(I,3)*ALINV(3,J)               
C        WRITE(6,702) I,J,FFF                                           
C 702    FORMAT(1H ,'ALTV(',I3,')*ALINV(',I3,')=',F20.15)               
  703    CONTINUE                                                       
C--*--GENERATE TRANSLATIONAL VECTOR                                     
      MM=0                                                              
      DO 500 I=-NEIBR,NEIBR                                             
        DO 510 J=-NEIBR,NEIBR                                           
          DO 520 K=-NEIBR,NEIBR                                         
            MM = MM + 1                                                 
            IRX(MM) = I                                                 
            IRY(MM) = I                                                 
            IRZ(MM) = I                                                 
                                      F1  = DFLOAT(I)                   
                                      F2  = DFLOAT(J)                   
                                      F3  = DFLOAT(K)                   
            RX(MM)  = ALTV(1,1)*F1 + ALTV(1,2)*F2 + ALTV(1,3)*F3        
            RY(MM)  = ALTV(2,1)*F1 + ALTV(2,2)*F2 + ALTV(2,3)*F3        
            RZ(MM)  = ALTV(3,1)*F1 + ALTV(3,2)*F2 + ALTV(3,3)*F3        
            RR(MM)  = SQRT( RX(MM)**2 + RY(MM)**2 + RZ(MM)**2 )         
  520     CONTINUE                                                      
  510   CONTINUE                                                        
  500 CONTINUE                                                          
      IF(MM.NE.NEIBRD) THEN                                             
        WRITE(6,*) ' MM.NE.NEIBRD MM,NEIBRD =',MM,NEIBRD                
      END IF                                                            
C     CALL HPSORT(NEIBRD,IRX,IRY,IRZ,RX,RY,RZ,RR) <------- STRESS       
      IF(IPRI.GE.2) THEN                                                
          WRITE(6,*) ' TRANSLATIONAL VECTOR'                            
          WRITE(6,600) (I,IRX(I),IRY(I),IRZ(I),                         
     &                        RX(I) ,RY(I) ,RZ(I) ,RR(I),I=1,NEIBRD)    
      END IF                                                            
  600 FORMAT((' ',4I5,4F12.6))                                          
C     WRITE (6,*) 'UNIVOL(IN LATTIC) = ',UNIVOL                         
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                
      SUBROUTINE GSTEP1(IPRE)                                           
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      KMAX=KNX+1                                                        
      KMAY=KNY+1                                                        
      KMAZ=KNZ+1                                                        
      IF (IPRE.NE.1) THEN                                               
        KMAX=KNX-1                                                      
        KMAY=KNY-1                                                      
        KMAZ=KNZ-1                                                      
      END IF                                                            
      KMAX2=2*KMAX+1                                                    
      KMAY2=2*KMAY+1                                                    
      KMAZ2=2*KMAZ+1                                                    
C     ## LARGEST WAVE VECTOR; CUT OFF WAVE NUMBER FOR WAVE FUNCTION ##  
      GMAX2=GMAX*2.0D0                                                  
C     //////////////////////////////////////////////                    
C     // GENERATION OF RECIPROCAL LATTICE VECTORS //                    
C     //                          STEP 1          //                    
C     //////////////////////////////////////////////                    
           B1X = RLTV(1,1)                                              
           B1Y = RLTV(2,1)                                              
           B1Z = RLTV(3,1)                                              
                           B2X = RLTV(1,2)                              
                           B2Y = RLTV(2,2)                              
                           B2Z = RLTV(3,2)                              
                                           B3X = RLTV(1,3)              
                                           B3Y = RLTV(2,3)              
                                           B3Z = RLTV(3,3)              
      MM=0                                                              
      DO 20 I1=1,KMAX2                                                  
                                      F1=FLOAT(I1-KMAX-1)               
      DO 20 I2=1,KMAY2                                                  
                                      F2=FLOAT(I2-KMAY-1)               
      DO 20 I3=1,KMAZ2                                                  
                                      F3=FLOAT(I3-KMAZ-1)               
         FX=  F1*B1X   +   F2*B2X   +   F3*B3X                          
         FY=  F1*B1Y   +   F2*B2Y   +   F3*B3Y                          
         FZ=  F1*B1Z   +   F2*B2Z   +   F3*B3Z                          
         FR=SQRT(FX*FX+FY*FY+FZ*FZ)                                     
             IF(FR.LE.GMAX2)  THEN                                      
             MM=MM+1                                                    
                 IG1(MM)=I1-KMAX-1                                      
                 IG2(MM)=I2-KMAY-1                                      
                 IG3(MM)=I3-KMAZ-1                                      
                              GX(MM)=FX                                 
                              GY(MM)=FY                                 
                              GZ(MM)=FZ                                 
                              GR(MM)=FR                                 
             END IF                                                     
20    CONTINUE                                                          
      WRITE(6,*) 'NUMBER OF GENERATED VECTORS=',MM                      
      IF(MM.NE.KG) WRITE(6,*) '**WARNING] KG MUST BE EQUAL TO  ',MM     
C-----------------------------------------------------------------------
      KG=MM                                                             
          IXMA=0                                                        
          IYMA=0                                                        
          IZMA=0                                                        
          DO 643 I=1,KG                                                 
          IXMA= MAX(IXMA,ABS(IG1(I)))                                   
          IYMA= MAX(IYMA,ABS(IG2(I)))                                   
          IZMA= MAX(IZMA,ABS(IG3(I)))                                   
643       CONTINUE                                                      
          IXMA = IXMA +1                                                
          IYMA = IYMA +1                                                
          IZMA = IZMA +1                                                
          WRITE(6,*) '**KX1 SHOULD BE =',IXMA                           
          WRITE(6,*) '**KY1 SHOULD BE =',IYMA                           
          WRITE(6,*) '**KZ1 SHOULD BE =',IZMA                           
C                                                                       
          KX1=IXMA                                                      
          KY1=IYMA                                                      
          KZ1=IZMA                                                      
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                
      SUBROUTINE GSTSCA(IPRE,EPP)                                       
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      REAL EPP(3,3),ERE(3,3),EPCI(3,3)
C     SXX, SXY, SXZ
C     SYX, SYY, SYZ
C     SZX, SZY, SZZ
C     X         Y         Z
C 1   ALTV(1,1),ALTV(2,1),ALTV(3,1)   AX1,AY1,AZ1
C 2   ALTV(1,2),ALTV(2,2),ALTV(3,2)   AX2,AY2,AZ2
C 3   ALTV(1,3),ALTV(2,3),ALTV(3,3)   AX3,AY3,AZ3
C
C     CALCULATE INVERSE MATRIX OF STRESS
      SXX=EPP(1,1)
      SXY=EPP(1,2)
      SXZ=EPP(1,3)
      SYY=EPP(2,2)
      SYZ=EPP(2,3)
      SZZ=EPP(3,3)
C
      ADJA=SXX*SYY*SZZ+2.0D0*SXY*SXZ*SYZ-SXZ*SXZ*SYY-SYZ*SYZ*SXX
     &    -SXY*SXY*SZZ
      ADJAI=1.0D0/ADJA
C
      ERE(1,1)=ADJAI*(SYY*SZZ-SYZ*SYZ)
      ERE(1,2)=ADJAI*(SYZ*SXZ-SZZ*SXY)
      ERE(1,3)=ADJAI*(SXY*SYZ-SXZ*SYY)
      ERE(2,1)=ERE(1,2)
      ERE(2,2)=ADJAI*(SZZ*SXX-SXZ*SXZ)
      ERE(2,3)=ADJAI*(SXZ*SXY-SXX*SYZ)
      ERE(3,1)=ERE(1,3)
      ERE(3,2)=ERE(2,3)
      ERE(3,3)=ADJAI*(SXX*SYY-SXY*SXY)
C     CHECK!
      DO 2000 I=1,3
      DO 2100 J=1,3
      EPCI(I,J)=EPP(I,1)*ERE(1,J)+EPP(I,2)*ERE(2,J)
     &         +EPP(I,3)*ERE(3,J)
 2100 CONTINUE
 2000 CONTINUE
C
      WRITE (6,*) 'CHECK! EPP*ERE'
      WRITE (6,*) EPCI(1,1),EPCI(1,2),EPCI(1,3)
      WRITE (6,*) EPCI(2,1),EPCI(2,2),EPCI(2,3)
      WRITE (6,*) EPCI(3,1),EPCI(3,2),EPCI(3,3)
C
      DO 1000 I=1,KNG
      GXX  =ERE(1,1)*GX(I)+ERE(1,2)*GY(I)+ERE(1,3)*GZ(I)
      GYY  =ERE(2,1)*GX(I)+ERE(2,2)*GY(I)+ERE(2,3)*GZ(I)
      GZZ  =ERE(3,1)*GX(I)+ERE(3,2)*GY(I)+ERE(3,3)*GZ(I)
      GX(I)=GXX
      GY(I)=GYY
      GZ(I)=GZZ 
      GR(I)=SQRT(GX(I)*GX(I)+GY(I)*GY(I)+GZ(I)*GZ(I))
 1000 CONTINUE
C
C     ERE(1,1)=1.0D0/EPP(1,1)                                           
C     ERE(2,2)=1.0D0/EPP(2,2)                                           
C     ERE(3,3)=1.0D0/EPP(3,3)                                           
      RETURN                                                            
      END                                                               
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7  
      SUBROUTINE SYMM(KOPR,NOPR,NBZTYP,NEIBRD,NFOUT,RX,RY,RZ)
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      DIMENSION NPPP(KNG,KO),NAAA(KATM,KO)
      DIMENSION OO(3,3,KO),TAA(3,KO)      
      DIMENSION RX(NEIBRD),RY(NEIBRD),RZ(NEIBRD)                        
!XOCL SUBPROCESSOR PT(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IQ=(PT,INDEX=1:KO,PART=BAND)
!XOCL LOCAL NPPP(:,/IQ),NAAA(:,/IQ),OO(:,:,/IQ),TAA(:,/IQ)
      EQUIVALENCE (NGPT,NPPP),(NAPT,NAAA),(OP,OO),(TAU,TAA)
C---------------------------------------------------------------------- 
      IF (NBZTYP.EQ.5) THEN                                             
          A = ABS(ALTV(1,1))*2.D0                                       
        ELSE IF(NBZTYP.EQ.6.OR.NBZTYP.EQ.7) THEN                        
          A = ABS(ALTV(1,1))                                            
          B = ABS(ALTV(2,2))                                            
        ELSE IF(NBZTYP.EQ.8 .OR. NBZTYP.EQ.10) THEN                     
          A = ABS(ALTV(1,1))                                            
          B = ABS(ALTV(3,3))                                            
        ELSE IF(NBZTYP.EQ.9.OR.NBZTYP.GE.11) THEN
          A = ABS(ALTV(1,1))                                            
          B = ABS(ALTV(2,2))                                            
          C = ABS(ALTV(3,3))                                            
      END IF                                                            
C---------------------------------------------                          
      CALL OPGR                                                         
     >     (NBZTYP,KOPR,A,B,C,NFOUT,                               
     <      NOPR)                                                
C--*--SYMMETRY PAIRS FOR G-POINTS ------------
      DO 100 I = 1,KG                                                   
        PX=GX(I)                                                        
        PY=GY(I)                                                        
        PZ=GZ(I)                                                        
!XOCL SPREAD DO /IQ
        DO 110 NO=1,NOPR                                                
          FX=OO(1,1,NO)*PX +OO(1,2,NO)*PY +OO(1,3,NO)*PZ                
          FY=OO(2,1,NO)*PX +OO(2,2,NO)*PY +OO(2,3,NO)*PZ                
          FZ=OO(3,1,NO)*PX +OO(3,2,NO)*PY +OO(3,3,NO)*PZ                
          DO 120 J = 1,KG                                               
            FF1 = ABS(GX(J)-FX)+ABS(GY(J)-FY)+ABS(GZ(J)-FZ)             
            IF (FF1.LE.1.D-5) THEN                                      
                NPPP(I,NO) = J                                          
                GOTO 110                                                
            END IF                                                      
  120     CONTINUE                                                      
          WRITE(NFOUT,130) I,NO                                         
          STOP                                                          
  110   CONTINUE                                                        
!XOCL END SPREAD
  100 CONTINUE                                                          
  130 FORMAT(' ','THERE IS NO PAIR FOR NG,NOP=',2I8)                    
C--*--SYMMETRY PAIRS FOR ATOMS                                          
      DDD = 1.0D-5                                                      
      DO 200 I = 1,KATM
        PX=CATX(I)                                                     
        PY=CATY(I)                                                     
        PZ=CATZ(I)                                                     
!XOCL SPREAD DO /IQ
        DO 210 NO=1,NOPR                                                
          FX=OO(1,1,NO)*PX+OO(1,2,NO)*PY+OO(1,3,NO)*PZ+TAA(1,NO)        
          FY=OO(2,1,NO)*PX+OO(2,2,NO)*PY+OO(2,3,NO)*PZ+TAA(2,NO)        
          FZ=OO(3,1,NO)*PX+OO(3,2,NO)*PY+OO(3,3,NO)*PZ+TAA(3,NO)        
          DO 220 J = 1,KATM
            DO 230 K = 1,NEIBRD                                         
              FFX = ABS( FX - CATX(J) - RX(K) )                        
              FFY = ABS( FY - CATY(J) - RY(K) )                        
              FFZ = ABS( FZ - CATZ(J) - RZ(K) )                        
              IF( (FFX.LE.DDD).AND.                                     
     &            (FFY.LE.DDD).AND.                                     
     &            (FFZ.LE.DDD)     ) THEN                               
                    NAAA(I,NO) = J                                      
                    GOTO 210                                            
              END IF                                                    
  230       CONTINUE                                                    
  220     CONTINUE                                                      
          WRITE(NFOUT,*) ' NO PAIR I, NO ',I,NO                         
          STOP                                                          
  210   CONTINUE                                                        
!XOCL END SPREAD
  200 CONTINUE                                                          
C--*--OUTPUT                                                            
      IF(IPRI.GE.2) THEN                                                
        WRITE(NFOUT,*) ' NGPT '                                         
        DO 310 I=1,20                                                   
          WRITE(NFOUT,*) ' NG =',I                                      
          WRITE(NFOUT,300) (NGPT(I,J),J=1,NOPR)                         
  310   CONTINUE                                                        
        DO 320 I=KG-20,KG                                               
          WRITE(NFOUT,*) ' NG =',I                                      
          WRITE(NFOUT,300) (NGPT(I,J),J=1,NOPR)                         
  320   CONTINUE                                                        
        WRITE(NFOUT,*) ' NAPT '                                         
        DO 330 I=1,KATM
          WRITE(NFOUT,*) ' NA =',I                                      
          WRITE(NFOUT,300) (NAPT(I,J),J=1,NOPR)                         
  330   CONTINUE                                                        
  300   FORMAT((8I8))                                                   
      END IF                                                            
      RETURN                                                            
      END                                                               
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7  
      SUBROUTINE SYMSCA(KOPR,NOPR,NBZTYP,NEIBRD,NFOUT,RX,RY,RZ,EPP)
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      DIMENSION TAA(3,KO)      
!XOCL SUBPROCESSOR PT(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IQ=(PT,INDEX=1:KO,PART=BAND)
!XOCL LOCAL TAA(:,/IQ)
      EQUIVALENCE (TAU,TAA)
      DIMENSION EPP(3,3)                                                 
C---------------------------------------------------------------------- 
!XOCL SPREAD DO /IQ
      DO 1000 I=1,NOPR                                                  
C     TAU(1,I)=EPP(1,1)*TAU(1,I)+EPP(1,2)*TAU(2,I)+EPP(1,3)*TAU(3,I)
C     TAU(2,I)=EPP(1,2)*TAU(1,I)+EPP(2,2)*TAU(2,I)+EPP(2,3)*TAU(3,I)
C     TAU(3,I)=EPP(1,3)*TAU(1,I)+EPP(2,3)*TAU(2,I)+EPP(3,3)*TAU(3,I)
      TAA(1,I)=EPP(1,1)*TAA(1,I)
      TAA(2,I)=EPP(2,2)*TAA(2,I)
      TAA(3,I)=EPP(3,3)*TAA(3,I)
 1000 CONTINUE                                                          
!XOCL END SPREAD
C     DO 1100 I=1,NOPR                                                  
C     WRITE (6,600) TAU(1,I),TAU(2,I),TAU(3,I)                          
C 600 FORMAT(1H ,'TAU(1,2,3) = ',3F12.6)                                
C1100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7  
      SUBROUTINE FORCES                                                 
     >          (KOPR,NOPR,FORC,FORCW)                                       
C FORCE SYMMETRIZATION                                                  
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      DIMENSION NAAA(KATM,KO)
      DIMENSION OO(3,3,KO)
!XOCL SUBPROCESSOR PT(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IQ=(PT,INDEX=1:KO,PART=BAND)
!XOCL LOCAL NAAA(:,/IQ),OO(:,:,/IQ)
      EQUIVALENCE (NAPT,NAAA),(OP,OO)
      DIMENSION FORC(KATM,3),FORCW(KATM,3)
C---------------------------------------------------------------------- 
      DENOM = 1.0D0/FLOAT(NOPR)                                         
      DO 100 IA = 1,KATM
        FORCW(IA,1) = 0.0D0                                             
        FORCW(IA,2) = 0.0D0                                             
        FORCW(IA,3) = 0.0D0                                             
  100 CONTINUE                                                          
!XOCL SPREAD DO /IQ
      DO 210 IOP = 1,NOPR                                             
        DO 200 IA= 1,KATM
          IAA= NAAA(IA,IOP)                                             
          F1 = FORC(IAA,1)                                              
          F2 = FORC(IAA,2)                                              
          F3 = FORC(IAA,3)                                              
          F4 = OO(1,1,IOP)*F1 +OO(2,1,IOP)*F2 +OO(3,1,IOP)*F3           
          F5 = OO(1,2,IOP)*F1 +OO(2,2,IOP)*F2 +OO(3,2,IOP)*F3           
          F6 = OO(1,3,IOP)*F1 +OO(2,3,IOP)*F2 +OO(3,3,IOP)*F3           
          FORCW(IA,1) = FORCW(IA,1) + F4                                
          FORCW(IA,2) = FORCW(IA,2) + F5                                
          FORCW(IA,3) = FORCW(IA,3) + F6                                
  200   CONTINUE                                                        
  210 CONTINUE                                                          
!XOCL END SPREAD SUM(FORCW)
      DO 300 IA = 1,KATM
        FORC(IA,1) = FORCW(IA,1) * DENOM                                
        FORC(IA,2) = FORCW(IA,2) * DENOM                                
        FORC(IA,3) = FORCW(IA,3) * DENOM                                
  300 CONTINUE                                                          
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                
      SUBROUTINE GSTEPF                                                 
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      DO 30 I=1,KG                                                      
      IF (IG1(I).LE.-1) IGF1(I)=IG1(I)+2*IFX-1                          
      IF (IG1(I).GE. 0) IGF1(I)=IG1(I)+1                                
      IF (IG2(I).LE.-1) IGF2(I)=IG2(I)+2*IFY-1                          
      IF (IG2(I).GE. 0) IGF2(I)=IG2(I)+1                                
      IF (IG3(I).LE.-1) IGF3(I)=IG3(I)+2*IFZ-1                          
      IF (IG3(I).GE. 0) IGF3(I)=IG3(I)+1                                
   30 CONTINUE                                                          
C     --------------------TETRAHEDRON VOLUME-------------------------   
      FFF =    RLTV(1,1)*(RLTV(2,2)*RLTV(3,3)-RLTV(3,2)*RLTV(2,3))      
     &      +  RLTV(2,1)*(RLTV(3,2)*RLTV(1,3)-RLTV(1,2)*RLTV(3,3))      
     &      +  RLTV(3,1)*(RLTV(1,2)*RLTV(2,3)-RLTV(2,2)*RLTV(1,3))      
      RVOL = ABS(FFF)                                                   
      DO 1000  K=1,KNZ2                                                 
      DO 1000  J=1,KNY2                                                 
      DO 1000  I=1,KNX2                                                 
      IGPO(I,J,K) = 0                                                   
 1000 CONTINUE                                                          
C     CALCULATION OF POINTER IGPO                                       
      DO 1100 I=1,KG                                                    
      IGPO(IG1(I)+KX1,IG2(I)+KY1,IG3(I)+KZ1)= I                         
 1100 CONTINUE                                                          
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                
      SUBROUTINE GSFSCA                                                 
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
C     --------------------TETRAHEDRON VOLUME-------------------------   
      FFF =    RLTV(1,1)*(RLTV(2,2)*RLTV(3,3)-RLTV(3,2)*RLTV(2,3))      
     &      +  RLTV(2,1)*(RLTV(3,2)*RLTV(1,3)-RLTV(1,2)*RLTV(3,3))      
     &      +  RLTV(3,1)*(RLTV(1,2)*RLTV(2,3)-RLTV(2,2)*RLTV(1,3))      
      RVOL = ABS(FFF)                                                   
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^          
      SUBROUTINE BASNUM(QWGT)                                           
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^          
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      DIMENSION QWGT(KNV3)
C     //////////////////////////////                                    
C     // CUT OFF FOR WAVEFUNCTION //                                    
C     //////////////////////////////                                    
      JGB=0                                                             
      JG2=0                                                             
      JJT=0                                                             
      GRVT=0.D0                                                         
      DO 1000 I=1,KV3                                                   
                      DX = VX(I)                                        
                      DY = VY(I)                                        
                      DZ = VZ(I)                                        
                      JJ = 0                                            
                      KK = 0                                            
                      GRVMX =  0.D0                                     
      DO 1100 J=1,KG                                                    
        GRVV = SQRT((DX+GX(J))**2+(DY+GY(J))**2+(DZ+GZ(J))**2)          
        IF(GRVV.LE.GMAX) THEN                                           
            JJ = JJ + 1                                                 
            NBASE(JJ,I)=J                                               
            IF(GRVMX.LT.GRVV) GRVMX = GRVV                              
        END IF                                                          
        IF(GRVV.LE.GMAX/2.0D0) THEN                                     
         KK = KK + 1                                                    
         NBMAT(KK,I)=J                                                  
        END IF                                                          
 1100 CONTINUE                                                          
      IBA(I) = JJ                                                       
      JJT    = JJT+JJ                                                   
      GRVT   = GRVT+GRVMX*QWGT(I)                                       
      IBA2(I)= KK                                                       
      WRITE(6,1200) I,JJ,KK,GRVT                                        
 1200 FORMAT(' ',' K = ',I4,' JJ = ',I5,' KK = ',I5,' GRV = ',F12.8)    
      IF(JJ.GT.JGB) JGB = JJ                                            
      IF(KK.GT.JG2) JG2 = KK                                            
 1000 CONTINUE                                                          
      WRITE(6,*) 'MAXIMUM NUMBER OF BASES FOR WAVE FUNCTION (KG1)=',JGB 
      WRITE(6,*) ' JJT = ',JJT,' MEAN GRV = ',GRVT                      
      WRITE(6,*) 'MAXIMUM NUMBER OF BASES FOR (GMAX/2)      (KG2)=',JG2 
C     <----- ATTN] KG1 WO SAITEIGI SHITE IRU]                           
      KG1=JGB                                                           
      NBMX = 0                                                          
      DO 2000 IK = 1,KV3                                                
        IIBA = IBA(IK)                                                  
        DO 2010 I= 1,IIBA                                               
          NBMX=MAX(NBASE(I,IK),NBMX)                                    
 2010   CONTINUE                                                        
 2000 CONTINUE                                                          
      WRITE(6,2020) NBMX                                                
 2020 FORMAT(' ','MAX NBMX  = ',I6)                                     
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                        
      SUBROUTINE FORM                                                   
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                        
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      DO 1000 IT = 1,KTYP                                               
        DO 1010 N = 1,KG                                                
          ZFM3(N,IT) = DCMPLX(0.0D0,0.0D0)                              
 1010   CONTINUE                                                        
 1000 CONTINUE                                                          
      DO 2000 IA=1,KATM                                                 
      DO 2010 I=1,KG                                                    
        FX = GX(I)                                                      
        FY = GY(I)                                                      
        FZ = GZ(I)                                                      
        ZFORT      =CDEXP(-ZI*(CATX(IA)*FX+CATY(IA)*FY+CATZ(IA)*FZ))    
*VOCL STMT,IF(90)
        IF(I.LE.KNBMX) ZFM2(I,IA) = ZFORT                               
        ZFM3(I,KFTYPE(IA)) = ZFM3(I,KFTYPE(IA)) + ZFORT                 
 2010 CONTINUE                                                          
 2000 CONTINUE                                                          
      RETURN                                                            
      END                                                               
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7  
      SUBROUTINE PCC(KTPCC,RHPC,                                      
     <                 SCHGPC,EPC,PCM)                
C     SUBROUTINE FOR PARTIAL CORE CORRECTION                            
C                         MODIFIED FROM      ATOMIC SEP.1990 MORIKAWA   
C                         REF.     S.G.LOUIE,S.FROYEN,AND M.L.COHEN,    
C                                  PHYS.REV.B26,1738(1982)              
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'
      INTEGER  I,IA,ITY,N,KTPCC(KTYP)                                   
      REAL     EPC,PCM,RS
      REAL     XH,H,H3,S,SCHGPC
      REAL     RHPC(KIZA,KTYP),CHGPC(KTYP),EPCP(KTYP),EMCP(KTYP)      
      DATA XH/32.0D0/                                                   
      H   = 1.D0/XH                                                     
      H3  = H/3.D0                                                      
                CO3=4.D0/3.D0                                           
C     WRITE (6,*) 'UNIVOL IN PCC =',UNIVOL                              
C/////////////  FOURIER TRANSFORMATION OF PARTIAL CHARGE DENSITY /////  
      IF (ITER.EQ.0) THEN                                               
      DO 1110 ITY=1,KTYP
      DO 1111 N=1,KIZA
      RHPC(N,ITY)=0.0D0
 1111 CONTINUE
 1110 CONTINUE
C----------------------------------------
        REWIND 27                                                       
        READ(27,1100) (I,RAD(N),N=1,KIZA)                               
 1100 FORMAT(3(I4,D22.14))                                              
        IF (I.NE.KIZA) THEN                                             
          WRITE (6,*) ' READ  IN PCC '                                  
          STOP                                                          
        END IF                                                          
        WRITE(6,*) I,' RAD(KIZA)= ',RAD(KIZA)                           
        DO 1000 ITY=1,KTYP                                              
          IF (KTPCC(ITY).EQ.0) GOTO 1000                                
            READ(27,*) CHGPC(ITY)                                       
            READ(27,1100) (I,RHPC(N,ITY),N=1,KIZA)                      
          IF (I.NE.KIZA) THEN                                           
            WRITE (6,*) ' READ  IN PCC  ITY= ',ITY                      
            STOP                                                        
          END IF                                                        
 1000 CONTINUE
      END IF
C--------------------------------
      DO 1001 ITY=1,KTYP
        WRITE(6,*) ITY,' RHPC(KIZA,ITY)= ',RHPC(KIZA,ITY)               
        S         = 0.D0                                                
        EPCP(ITY) = 0.D0                                                
        EMCP(ITY) = 0.D0                                                
        DO 1150 N=1,KIZA                                                
          IF (RHPC(N,ITY).GE.1.D-40) THEN                               
            RS = (3.D0/RHPC(N,ITY))**(1.D0/3.D0)*RAD(N)                 
            EPCP(ITY) = EPCP(ITY)-(0.458D0/RS+0.44D0/(RS+7.8D0))        
     &                                              *RHPC(N,ITY)*OMO(N) 
            EMCP(ITY) = EMCP(ITY)-(0.458D0*CO3/RS+(0.44D0*CO3*RS +      
     &                       3.432D0)/(RS+7.8D0)**2)*RHPC(N,ITY)*OMO(N) 
          END IF                                                        
          S           = S + RHPC(N,ITY)*OMO(N)                          
 1150   CONTINUE                                                        
        S         =     S    *H3                                        
        EPCP(ITY) = EPCP(ITY)*H3                                        
        EMCP(ITY) = EMCP(ITY)*H3                                        
        WRITE (6,*) ' CHGPC(ITY)  = ',CHGPC(ITY)                        
        IF (ABS(S-CHGPC(ITY)).GE.1.D-7) THEN                            
          WRITE (6,*) ITY,'PARTIAL CORE CHARGE  S-CHGPC='               
     /                                                    ,S-CHGPC(ITY) 
          STOP                                                          
        END IF                                                          
C------------------------------------
        DO 1300 I=1,KG
        RHPCG(I,ITY)=0.0D0
        RRPCG(I,ITY)=0.0D0
 1300   CONTINUE
        DO 1200 N=1,KIZA
        R=RAD(N)
        FAC1=OMO(N)*RHPC(N,ITY)*H3/UNIVOL
        FAC2=OMO(N)*RHPC(N,ITY)*R*H3/UNIVOL
        DO 1210 I=1,KG
        XX(I)=R*GR(I)
 1210   CONTINUE
        CALL DNV(0,KG,XX,Y11)
        CALL DNV(1,KG,XX,Y22)
        DO 1220 I=1,KG
        RHPCG(I,ITY)=RHPCG(I,ITY)+FAC1*Y11(I)
        RRPCG(I,ITY)=RRPCG(I,ITY)-FAC2*Y22(I)
 1220   CONTINUE
 1200   CONTINUE
 1001 CONTINUE                                                          
      SCHGPC = 0.D0                                                     
      EPC = 0.D0                                                        
      PCM = 0.0D0                                                       
      DO 1410 IA = 1,KATM                                               
        IF (KTPCC(KFTYPE(IA)).EQ.0) GOTO 1410                           
        SCHGPC = SCHGPC+CHGPC(KFTYPE(IA))                               
        EPC    = EPC   +EPCP (KFTYPE(IA))                               
        PCM    = PCM   +EMCP (KFTYPE(IA))                               
 1410 CONTINUE                                                          
      WRITE (6,*) 'NO-OVERLAP P.C.(USED) EPC = ',EPC                    
      WRITE (6,*) 'NO-OVERLAP P.C.(USED) PCM = ',PCM                    
      WRITE (6,*) 'TOTAL PARTIAL CHARGE    SCHGPC = ',SCHGPC            
C///////////////////////////////////////////////////////////////////    
        ZRHPC(1)=DCMPLX(0.D0,0.D0)                                      
        ZRRPC(1)=DCMPLX(0.D0,0.D0)                                      
        RG(1)=1.0D0                                                     
      DO 1500 I=2,KG                                                    
        ZRHPC(I)=DCMPLX(0.D0,0.D0)                                      
        ZRRPC(I)=DCMPLX(0.D0,0.D0)                                      
        RG(I)=1.0D0/GR(I)                                               
 1500 CONTINUE                                                          
      DO 2000 IT=1,KTYP                                                 
        IF (KTPCC(IT).EQ.0) GOTO 2000                                   
        DO 2100 I=1,KG                                                  
          ZRHPC(I) = ZRHPC(I) +  ZFM3(I,IT)*RHPCG(I,IT)                 
          ZRRPC(I) = ZRRPC(I) +  ZFM3(I,IT)*RRPCG(I,IT)*RG(I)           
 2100   CONTINUE                                                        
 2000 CONTINUE                                                          
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                    
      SUBROUTINE INTCHG(ALFA)                                           
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                    
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      DO 128 I=1,KG                                                     
      ZCHG(I) = DCMPLX(0.0D0,0.0D0)                                     
  128 CONTINUE                                                          
              DO 15 IT=1,KTYP                                           
              DO 10 I=1,KG                                              
              ZCHG(I) = ZCHG(I) +  ZFM3(I,IT)*                          
     &                  ACHG(IT)*EXP(-GR(I)*GR(I)*0.25D0/ALFA)/UNIVOL   
10    CONTINUE                                                          
15    CONTINUE                                                          
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                         
      SUBROUTINE FORLOC(KTPCC)                             
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      DIMENSION ZFPC(6)
      INTEGER   KTPCC(KTYP)                                             
C     FORCE2 AND PRESSURE ADDITION SSN(KNG1)                            
      DO 161 IA=1,KATM                                                  
      DO 162 I=1,KG                                                     
        FX = GX(I)                                                      
        FY = GY(I)                                                      
        FZ = GZ(I)                                                      
        ZFORT     =CDEXP(-ZI*(CATX(IA)*FX+CATY(IA)*FY+CATZ(IA)*FZ))    
        ZFTMP     =ZFORT  *PSC(I,KFTYPE(IA))*DCONJG(ZCHG(I))        
      ZFORC2(IA,1)=ZFORC2(IA,1)+GX(I)*ZFTMP                             
      ZFORC2(IA,2)=ZFORC2(IA,2)+GY(I)*ZFTMP                             
      ZFORC2(IA,3)=ZFORC2(IA,3)+GZ(I)*ZFTMP                             
  162 CONTINUE                                                          
      IF (KTPCC(KFTYPE(IA)).EQ.1) THEN                                  
          DO 1000 N=1,6                                                 
            ZFPC(N)=DCMPLX(0.0D0,0.0D0)                                 
 1000     CONTINUE                                                      
          DO 200 I=1,KG                                                 
C     4/20, 1999, MODIFIED SERIOUS BUG
            FX = GX(I)                                                      
            FY = GY(I)                                                      
            FZ = GZ(I)                                                      
            ZFORT =CDEXP(-ZI*(CATX(IA)*FX+CATY(IA)*FY+CATZ(IA)*FZ))    
            ZFTMP =ZFORT  *RHPCG(I,KFTYPE(IA))*DCONJG(ZVXC(I))      
            ZFPC(1)=ZFPC(1)+GX(I)*ZFTMP                                 
            ZFPC(2)=ZFPC(2)+GY(I)*ZFTMP                                 
            ZFPC(3)=ZFPC(3)+GZ(I)*ZFTMP                                 
            ZFTMP =ZFORT  *RHPCG(I,KFTYPE(IA))*DCONJG(ZVXCPC(I))    
            ZFPC(4)=ZFPC(4)+GX(I)*ZFTMP                                 
            ZFPC(5)=ZFPC(5)+GY(I)*ZFTMP                                 
            ZFPC(6)=ZFPC(6)+GZ(I)*ZFTMP                                 
  200     CONTINUE                                                      
          ZFORC2(IA,1)=ZFORC2(IA,1)+ZFPC(1)-ZFPC(4)                     
          ZFORC2(IA,2)=ZFORC2(IA,2)+ZFPC(2)-ZFPC(5)                     
          ZFORC2(IA,3)=ZFORC2(IA,3)+ZFPC(3)-ZFPC(6)                     
C          IF(MOD(ITER,50).EQ.0) THEN                                   
C              WRITE(6,1100)  IA,ZFPC(1),ZFPC(2),ZFPC(3),               
C     &                          ZFPC(4),ZFPC(5),ZFPC(6)                
C 1100         FORMAT(I4,6D12.4,/,4X,6D12.4)                            
C          END IF                                                       
      END IF                                                            
C     COMBINE THE FORCE1 AND FORCE2.                                    
  161 CONTINUE                                                          
      DO 261 IA=1,KATM
      ZFORC2(IA,1)=UNIVOL*ZI*ZFORC2(IA,1) + DCMPLX(FFF1(IA,1),0.0D0)    
      ZFORC2(IA,2)=UNIVOL*ZI*ZFORC2(IA,2) + DCMPLX(FFF1(IA,2),0.0D0)    
      ZFORC2(IA,3)=UNIVOL*ZI*ZFORC2(IA,3) + DCMPLX(FFF1(IA,3),0.0D0)    
  261 CONTINUE
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                       
      SUBROUTINE CONV2(ETONEW,ETOOLD,IMD,ITERMD,CMIX)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                       
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      IF (ITER.EQ.IMD+1) WRITE (6,*) 'MIXING RATE = ',CMIX
      IF (ITER.GE.IMD+1) PMIX(ITER)=CMIX
      ZFFF=DCMPLX(0.0D0,0.0D0)                                          
      PPMIX = PMIX(ITER)                                                
      QQMIX = 1.D0 - PPMIX                                              
        DO 1100 I=1,KG                                                  
      ZUU1=ZCHG(I)-ZCHGO(I)                                             
      ZFFF=ZFFF+DCONJG(ZUU1)*ZUU1                                       
        ZCHO1 = ZCHGO(I)                                             
        ZCHGO(I) = ZCHG(I)                                              
        ZCHG(I)  = ZCHO1*PPMIX + ZCHGO(I)*QQMIX                      
 1100 CONTINUE                                                          
      DCHG=DSQRT(DREAL(ZFFF))                                           
 5000 CONTINUE                                                          
      IF(ITER.NE.1) THEN                                                
                    EPP2=ABS(ETONEW-ETOOLD)                             
                                IF(ETOOLD.LT.ETONEW) THEN               
C-------------------------------WRITE(6,*) '**********************'     
                                WRITE(6,*) '>> ETOOLD.LT.ETONEW <<'     
C-------------------------------WRITE(6,*) '**********************'     
                                WRITE(6,*) 'ETONEW-ETTOOLD=',           
     &                                      ETONEW-ETOOLD               
                                EPP2=1.0D0                              
                                END IF                                  
                    ETOOLD=ETONEW                                       
        ELSE                                                            
              EPP2  =1.0D5                                              
              ETOOLD=ETONEW                                             
        END IF                                                          
        EDEV=27200.0D0*EPP2                                             
        WRITE(6,610)  ITER,EPP2,EDEV,DCHG                               
  610   FORMAT(1H ,' ITER=',I3,' ET(H)=',D15.7,' ET(M)=',F10.6,' DC='   
     &            ,D15.7)                                               
      IF( ITER.LE.ITEMAX  ) THEN                                        
C         IF(EPP2.LE.CONV) GOTO 2000                                    
          IF(ITER.NE.ITEMAX) THEN                                       
             IFLAG=1                                                    
             GO TO 9999                                                 
          END IF                                                        
      END IF                                                            
      WRITE(6,*) '////////////////////////////////////////////'         
      WRITE(6,*) '// UNLUCKY! CONVERGENCE WAS NOT ACHIEVED! //'         
      WRITE(6,*) '////////////////////////////////////////////'         
      GO TO 2001                                                        
 2000 CONTINUE                                                          
      WRITE(6,*) '///////////////////////////////////////////////'      
      WRITE(6,*) '// CONGRATULATION! CONVERGENCE WAS ACHIEVED! //'      
      WRITE(6,*) '///////////////////////////////////////////////'      
 2001 CONTINUE                                                          
      IFLAG=0                                                           
 9999 CONTINUE                                                          
      RETURN                                                            
      END                                                               
C================================================================       
      SUBROUTINE EVIN(IREC8,IREC9,EPC)
      IMPLICIT REAL(A-H,O-Y)
      IMPLICIT COMPLEX(Z)
      INCLUDE 'PACVPP'
      DIMENSION DVS(KIZA),DVP(KIZA),DVD(KIZA)
      REWIND 28
      REWIND 29
      REWIND 30
      REWIND 31
      REWIND 32
      REWIND 33
      REWIND 34
      REWIND 35
      REWIND 36
      REWIND 37
      READ (28,*) ITER,EPC,WS,WP,WD
      READ (29) SNL
      READ (30) SNL2
      READ (31) SNL3
      READ (32) ZAJ
      READ (33) ZFBB,ZFBB2
      READ (34) EKO
      READ (35) ZCHG
      READ (36) VEL
      READ (37) ZRHPC
      WRITE (6,*) 'EVIN  ITER = ',ITER,EPC
      DO 1000 I=NBD1,NBD2
      DO 1010 K=1,KNV3
      WRITE (6,600) EKO(I,K)
  600 FORMAT(1H ,4F12.6)
 1010 CONTINUE
 1000 CONTINUE
C     NSEK = 2 BECAUSE DVS,DVP INCLUDES RAD(N).        
      NSEK=2                                           
      REWIND 15                                                         
      DO 100 ITY = 1,KTYP                                               
        IF (NLSPD(ITY).EQ.1) THEN
          READ (15,*) MESHR,NMES,DX,WS(ITY),WP(ITY)
     &             ,RAD,DVS,DVP
          WD(ITY)=0.0D0
          DO 111 N=1,MESHR                                                
            WVS(N,ITY)=DVS(N)                                             
            WVP(N,ITY)=DVP(N)                                             
            WVD(N,ITY)=0.0D0
  111 CONTINUE
        ELSE
          READ (15,*) MESHR,NMES,DX,WS(ITY),WP(ITY),WD(ITY)
     &               ,RAD,DVS,DVP,DVD
          DO 110 N=1,MESHR                                                
            WVS(N,ITY)=DVS(N)                                             
            WVP(N,ITY)=DVP(N)                                             
            WVD(N,ITY)=DVD(N)
  110 CONTINUE                                                        
        END IF
  100 CONTINUE                                                          
      RETURN
      END
C================================================================       
      SUBROUTINE EVOUT(IREC8,IREC9,ALFA,EPC)
      IMPLICIT REAL(A-H,O-Y)
      IMPLICIT COMPLEX(Z)
      INCLUDE 'PACVPP'
      REWIND 11
      REWIND 20
      REWIND 28
      REWIND 29
      REWIND 30
      REWIND 31
      REWIND 32
      REWIND 33
      REWIND 34
      REWIND 35
      REWIND 36
      REWIND 37
      WRITE (20,300) EF
  300 FORMAT(1H ,'E-K CURVE EF =',F12.6)
      DO 100 NNN=1,KV3
      WRITE(20,303) VX(NNN),VY(NNN),VZ(NNN)
  303 FORMAT((1H ,3(F10.6,2X)))
      WRITE(20,302) (EKO(I,NNN),I=1,KEG)
  302 FORMAT((1H ,5(F10.6,2X)))
  100 CONTINUE
C
      WRITE (28,*) ITER,EPC,WS,WP,WD
      WRITE (29) SNL
      WRITE (30) SNL2
      WRITE (31) SNL3
      WRITE (32) ZAJ
      WRITE (33) ZFBB,ZFBB2
      WRITE (34) EKO
      WRITE (35) ZCHG
      WRITE (36) VEL
      WRITE (37) ZRHPC
C'''''WRITE  11!'''''''''''''''''''''''''''''''''''''''''''''''''''
      WRITE(11,399) ITEMAX,PINIT,CONV,GMAX,ICONT
  399 FORMAT(' ',I6,F8.4,D12.4,F8.4,I4)
      DO 400 I=1,3
        WRITE(11,500) ALTV(1,I),ALTV(2,I),ALTV(3,I)
  400 CONTINUE
      WRITE(11,*) KCOTYP,' COORDINATES  0:NORMALIZED  1:CARTESIAN '
      IF(KCOTYP.EQ.0) THEN
          DO 410 IA=1,KATM
            PPOS(IA)  =   ALINV(1,1)*CATX(IA)
     &                  + ALINV(1,2)*CATY(IA) + ALINV(1,3)*CATZ(IA)
            QPOS(IA)  =   ALINV(2,1)*CATX(IA)
     &                  + ALINV(2,2)*CATY(IA) + ALINV(2,3)*CATZ(IA)
            RPOS(IA)  =   ALINV(3,1)*CATX(IA)
     &                  + ALINV(3,2)*CATY(IA) + ALINV(3,3)*CATZ(IA)
  410     CONTINUE
          WRITE(11,500)(PPOS(IA),QPOS(IA),RPOS(IA),IA=1,KATM)
        ELSE
          WRITE(11,500)(CATX(IA),CATY(IA),CATZ(IA),IA=1,KATM)
      END IF
      NNATM=0
      DO 420 IT=1,KTYP
        NNATM=NNATM+IATOM(IT)
        WRITE(11,*) IATOM(IT),NLSPD(IT)
        WRITE(11,502) AICHG(NNATM),ALFA
        WRITE(11,502) ACHG(IT),AC(IT,1),AC(IT,2),BC(IT,1),BC(IT,2)
  420 CONTINUE
  500 FORMAT(3(F20.10))
  502 FORMAT(5(F15.8))
C'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
      RETURN
      END
C----------------------------------------------------------------
      SUBROUTINE EVOU2(ALFA)
      IMPLICIT REAL(A-H,O-Y)
      IMPLICIT COMPLEX(Z)
      INCLUDE 'PACVPP'
      REWIND 11
      REWIND 20
      WRITE (20,300) EF
  300 FORMAT(1H ,'E-K CURVE EF =',F12.6)
      DO 100 NNN=1,KV3
      WRITE(20,303) VX(NNN),VY(NNN),VZ(NNN)
  303 FORMAT((1H ,3(F10.6,2X)))
      WRITE(20,302) (EKO(I,NNN),I=1,KEG)
  302 FORMAT((1H ,5(F10.6,2X)))
  100 CONTINUE
C'''''WRITE  11!'''''''''''''''''''''''''''''''''''''''''''''''''''
      WRITE(11,399) ITEMAX,PINIT,CONV,GMAX,ICONT
  399 FORMAT(' ',I6,F8.4,D12.4,F8.4,I4)
      DO 400 I=1,3
        WRITE(11,500) ALTV(1,I),ALTV(2,I),ALTV(3,I)
  400 CONTINUE
      WRITE(11,*) KCOTYP,' COORDINATES  0:NORMALIZED  1:CARTESIAN '
      IF(KCOTYP.EQ.0) THEN
          DO 410 IA=1,KATM
            PPOS(IA)  =   ALINV(1,1)*CATX(IA)
     &                  + ALINV(1,2)*CATY(IA) + ALINV(1,3)*CATZ(IA)
            QPOS(IA)  =   ALINV(2,1)*CATX(IA)
     &                  + ALINV(2,2)*CATY(IA) + ALINV(2,3)*CATZ(IA)
            RPOS(IA)  =   ALINV(3,1)*CATX(IA)
     &                  + ALINV(3,2)*CATY(IA) + ALINV(3,3)*CATZ(IA)
  410     CONTINUE
          WRITE(11,500)(PPOS(IA),QPOS(IA),RPOS(IA),IA=1,KATM)
        ELSE
          WRITE(11,500)(CATX(IA),CATY(IA),CATZ(IA),IA=1,KATM)
      END IF
      NNATM=0
      DO 420 IT=1,KTYP
        NNATM=NNATM+IATOM(IT)
        WRITE(11,*) IATOM(IT),NLSPD(IT)
        WRITE(11,502) AICHG(NNATM),ALFA
        WRITE(11,502) ACHG(IT),AC(IT,1),AC(IT,2),BC(IT,1),BC(IT,2)
  420 CONTINUE
  500 FORMAT(3(F20.10))
  502 FORMAT(5(F15.8))
C'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
      RETURN
      END
C------------------------------------------------------------           
      SUBROUTINE XCFFT(III,SCHGPC)                         
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
C-----ARRAYS FOR MFFT--------------------------------------------       
      DIMENSION IWL(8*IFX2+28),IWM(8*IFY2+28),IWN(8*IFZ2+28)            
     &         ,IWORK(2*IFX2)                                           
      DIMENSION ZRB(IFX2,IFY2,IFZ2)
C============================================================           
      EQUIVALENCE (ZV1D(1),ZRB(1,1,1))                                   
C============================================================           
      KFT1 =  IFX2                                                      
      KFT2 =  IFY2                                                      
      KFT3 =  IFZ2                                                      
      KSUM=KFT1*KFT2*KFT3                                               
      KVOL=(KFT1-1)*KFT2*KFT3                                           
                THIRD=1.D0/3.D0                                         
                ALFF=0.7D0                                              
                CO1=1.5D0*ALFF*(3.D0/PAI)**THIRD                        
                CO2=(0.75D0/PAI)**THIRD                                 
                CO3=4.D0/3.D0                                           
      VCEL=ABS(UNIVOL/FLOAT(KVOL))                                      
C*****---- SET UP C3FFT (FAST FOURIER TRANSFORMATION) -----             
      CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                  
     &          ,0,0,1,IWORK,IERR)                                      
      IF (IERR.NE.0) THEN                                               
          WRITE (6,*) ' C3FFT(SET UP)]  IERR = ',IERR                   
          STOP                                                          
      END IF                                                            
      DO 111 I=1,KSUM                                                   
      ZV1D(I)=DCMPLX(0.0D0,0.0D0)                                        
  111 CONTINUE                                                          
      IF (III.EQ.3) THEN                                                
          DO 110 I=1,KG                                                 
            L1=IGF1(I)                                                  
            L2=IGF2(I)                                                  
            L3=IGF3(I)                                                  
            ZRB(L1,L2,L3)=        ZRHPC(I)                              
  110     CONTINUE                                                      
        ELSE                                                            
          DO 210 I=1,KG                                                 
            L1=IGF1(I)                                                  
            L2=IGF2(I)                                                  
            L3=IGF3(I)                                                  
            ZRB(L1,L2,L3)=ZCHG(I)+ZRHPC(I)                              
  210     CONTINUE                                                      
      END IF                                                            
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----                    
      CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                  
     &          ,0,1,1,IWORK,IERRI)                                     
      IF (IERRI.NE.0) THEN                                              
          WRITE (6,*) ' C3FFT(INVERSE)]  IERR = ',IERRI                 
          STOP                                                          
      END IF                                                            
      DO 1000 I=1,KSUM                                                  
        CHGB1(I) = DREAL(ZV1D(I))                                          
 1000 CONTINUE                                                          
      DO 20 I=1,KSUM                                                    
        IF( CHGB1(I).LE.0.0D0 ) THEN                                      
          IF ( MOD(I,IFX2).NE.0 ) THEN                                  
C              IF ( CHGB1(I).LE.-1.0D-5 ) THEN                            
C                WRITE(6,610) I,CHGB1(I)                                  
C  610           FORMAT(1H ,'**WARNING CHG.DEN<0.0 AT ',                
C     &                                  I5,2D15.7,'***')               
C              END IF                                                   
              CHGB1(I) = 1.D-40                                           
            ELSE                                                        
              CHGB1(I) = 1.D-40                                           
          END IF                                                        
        END IF                                                          
   20 CONTINUE                                                          
      S = 0.D0                                                          
      DO 1010 I=1,KSUM                                                  
        S = S+CHGB1(I)                                                    
 1010 CONTINUE                                                          
      S=S*VCEL - SCHGPC                                                 
      IF (MOD(ITER,10).EQ.0) THEN                                       
          WRITE (6,*) 'REAL TOTAL CHARGE = ',S,' IN XCFFT'              
      END IF                                                            
C*****<<< WIGNER INTERPOLATION >>>                                      
      IF((III.EQ.1).OR.(III.EQ.3)) THEN                                 
          DO 1020 I=1,KSUM                                              
            RS=CO2*( 1.0D0/CHGB1(I) )**THIRD                              
            RHO= -0.458D0*CO3/RS-(0.44D0*CO3*RS + 3.432D0)/             
     &                                          (RS+7.8D0)**2           
            ZV1D(I)=DCMPLX(RHO,0.0D0)                                    
 1020     CONTINUE                                                      
        ELSE                                                            
          DO 1030 I=1,KSUM                                              
            RS=CO2*( 1.0D0/CHGB1(I) )**THIRD                              
            RHO= -0.458D0/RS-0.44D0/(RS+7.8D0)                          
            ZV1D(I)=DCMPLX(RHO,0.0D0)                                    
 1030     CONTINUE                                                      
      END IF                                                            
C*****---- FAST FOURIER TRANSFORMATION -----                            
      CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                  
     &          ,0,-1,1,IWORK,IERR)                                     
      IF (IERR.NE.0) THEN                                               
          WRITE (6,*) ' C3FFT(FFT)]  IERR = ',IERR                      
          STOP                                                          
      END IF                                                            
      DENOM=1.0D0/DBLE(KVOL)                                            
C----------------------------------------------------
      IF (III.EQ.1) THEN                                                
        DO 120 I=1,KG                                                      
        L1=IGF1(I)                                                      
        L2=IGF2(I)                                                      
        L3=IGF3(I)                                                      
        ZVXC(I)=ZRB(L1,L2,L3)*DENOM                                      
  120 CONTINUE                                                          
        ELSE IF (III.EQ.2) THEN                                         
        DO 121 I=1,KG                                                      
        L1=IGF1(I)                                                      
        L2=IGF2(I)                                                      
        L3=IGF3(I)                                                      
        ZEXC(I)=ZRB(L1,L2,L3)*DENOM                                      
  121 CONTINUE                                                          
        ELSE IF (III.EQ.3) THEN                                         
        DO 122 I=1,KG                                                      
        L1=IGF1(I)                                                      
        L2=IGF2(I)                                                      
        L3=IGF3(I)                                                      
        ZVXCPC(I)=ZRB(L1,L2,L3)*DENOM                                      
  122 CONTINUE                                                          
      END IF                                                            
      RETURN                                                            
      END                                                               
C------------------------------------------------------------           
      SUBROUTINE XSTPC(III,SCHGPC)           
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'
C  STRESS FOR PARTIAL CORE CORRECTION
C-----ARRAYS FOR MFFT--------------------------------------------       
      DIMENSION IWL(8*IFX2+28),IWM(8*IFY2+28),IWN(8*IFZ2+28)            
     &         ,IWORK(2*IFX2)                                           
      DIMENSION ZRB(IFX2,IFY2,IFZ2),ZRX(KFFT)       
     &         ,ZRY(KFFT),ZRZ(KFFT)
C============================================================           
      EQUIVALENCE (ZV1D(1),ZRB(1,1,1))                                   
C============================================================           
      ICOUNT=0                                                          
      KFT1 =  IFX2                                                      
      KFT2 =  IFY2                                                      
      KFT3 =  IFZ2                                                      
      KSUM=KFT1*KFT2*KFT3                                               
      KVOL=(KFT1-1)*KFT2*KFT3                                           
                THIRD=1.D0/3.D0                                         
                ALFF=0.7D0                                              
                CO1=1.5D0*ALFF*(3.D0/PAI)**THIRD                        
                CO2=(0.75D0/PAI)**THIRD                                 
                CO3=4.D0/3.D0                                           
      VCEL=ABS(UNIVOL/FLOAT(KVOL))                                      
C*****---- SET UP C3FFT (FAST FOURIER TRANSFORMATION) -----             
      CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                  
     &          ,0,0,1,IWORK,IERR)                                      
      IF (IERR.NE.0) THEN                                               
          WRITE (6,*) ' C3FFT(SET UP)]  IERR = ',IERR                   
          STOP                                                          
      END IF                                                            
C----- X ---------------------------                                    
      DO 111 I=1,KSUM                                                   
      ZV1D(I)=DCMPLX(0.0D0,0.0D0)                                        
  111 CONTINUE                                                          
          DO 110 I=1,KG                                                 
            L1=IGF1(I)                                                  
            L2=IGF2(I)                                                  
            L3=IGF3(I)                                                  
            ZRB(L1,L2,L3)=        ZRRPC(I)*GX(I)*GX(I)                  
  110     CONTINUE                                                      
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----                    
      CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                  
     &          ,0,1,1,IWORK,IERRI)                                     
      IF (IERRI.NE.0) THEN                                              
          WRITE (6,*) ' C3FFT(INVERSE)]  IERR = ',IERRI                 
          STOP                                                          
      END IF                                                            
      DO 1005 I=1,KSUM                                                  
        ZRX(I) = ZV1D(I)                                                 
 1005 CONTINUE                                                          
C----- Y ---------------------------                                    
      DO 113 I=1,KSUM                                                   
      ZV1D(I)=DCMPLX(0.0D0,0.0D0)                                        
  113 CONTINUE                                                          
          DO 134 I=1,KG                                                 
            L1=IGF1(I)                                                  
            L2=IGF2(I)                                                  
            L3=IGF3(I)                                                  
            ZRB(L1,L2,L3)=        ZRRPC(I)*GY(I)*GY(I)                  
  134     CONTINUE                                                      
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----                    
      CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                  
     &          ,0,1,1,IWORK,IERRI)                                     
      IF (IERRI.NE.0) THEN                                              
          WRITE (6,*) ' C3FFT(INVERSE)]  IERR = ',IERRI                 
          STOP                                                          
      END IF                                                            
      DO 1006 I=1,KSUM                                                  
        ZRY(I) = ZV1D(I)                                                 
 1006 CONTINUE                                                          
C----- Z ---------------------------                                    
      DO 114 I=1,KSUM                                                   
      ZV1D(I)=DCMPLX(0.0D0,0.0D0)                                        
  114 CONTINUE                                                          
          DO 124 I=1,KG                                                 
            L1=IGF1(I)                                                  
            L2=IGF2(I)                                                  
            L3=IGF3(I)                                                  
            ZRB(L1,L2,L3)=        ZRRPC(I)*GZ(I)*GZ(I)                  
  124     CONTINUE                                                      
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----                    
      CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                  
     &          ,0,1,1,IWORK,IERRI)                                     
      IF (IERRI.NE.0) THEN                                              
          WRITE (6,*) ' C3FFT(INVERSE)]  IERR = ',IERRI                 
          STOP                                                          
      END IF                                                            
      DO 1007 I=1,KSUM                                                  
        ZRZ(I) = ZV1D(I)                                                 
 1007 CONTINUE                                                          
C-------------------------------                                        
      DO 112 I=1,KSUM                                                   
      ZV1D(I)=DCMPLX(0.0D0,0.0D0)                                        
  112 CONTINUE                                                          
          DO 210 I=1,KG                                                 
            L1=IGF1(I)                                                  
            L2=IGF2(I)                                                  
            L3=IGF3(I)                                                  
            ZRB(L1,L2,L3)=ZCHG(I)+ZRHPC(I)                             
  210     CONTINUE                                                      
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----                    
      CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                  
     &          ,0,1,1,IWORK,IERRI)                                     
      IF (IERRI.NE.0) THEN                                              
          WRITE (6,*) ' C3FFT(INVERSE)]  IERR = ',IERRI                 
          STOP                                                          
      END IF                                                            
      DO 1000 I=1,KSUM                                                  
        CHGB1(I) = DREAL(ZV1D(I))                                          
 1000 CONTINUE                                                          
      DO 20 I=1,KSUM                                                    
        IF( CHGB1(I).LE.0.0D0 ) THEN                                      
          IF ( MOD(I,IFX2).NE.0 ) THEN                                  
              IF ( CHGB1(I).LE.-1.0D-5 ) THEN                             
                ICOUNT=ICOUNT+1                                        
C               WRITE(6,610) I,CHGB1(I)                                   
C 610           FORMAT(1H ,'**WARNING CHG.DEN<0.0 AT ',                 
C    &                                  I5,2D15.7,'***')                
              END IF                                                    
              CHGB1(I) = 1.D-40                                           
            ELSE                                                        
              CHGB1(I) = 1.D-40                                           
          END IF                                                        
        END IF                                                          
   20 CONTINUE                                                          
      S = 0.D0                                                          
      DO 1010 I=1,KSUM                                                  
        S = S+CHGB1(I)                                                    
 1010 CONTINUE                                                          
      S=S*VCEL - SCHGPC                                                 
          WRITE (6,*) 'REAL TOTAL CHARGE = ',S,' IN XSTPC'              
C  -- X ------------                                                    
          DO 1020 I=1,KSUM                                              
            RS=CO2*( 1.0D0/CHGB1(I) )**THIRD                              
            RH1   = -0.458D0*CO3/RS-(0.44D0*CO3*RS + 3.432D0)/          
     &                                          (RS+7.8D0)**2           
            RH2   = -0.458D0/RS-0.44D0/(RS+7.8D0)                       
C                                                                       
            RHV(I)=RH1 - RH2 
            ZV1D(I)=ZRX(I)*(RHV(I))/CHGB1(I)
 1020     CONTINUE                                                      
C*****---- FAST FOURIER TRANSFORMATION -----                            
      CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                  
     &          ,0,-1,1,IWORK,IERR)                                     
      IF (IERR.NE.0) THEN                                               
          WRITE (6,*) ' C3FFT(FFT)]  IERR = ',IERR                      
          STOP                                                          
      END IF                                                            
      DENOM=1.0D0/DBLE(KVOL)                                            
      DO 30 I=1,KG                                                      
        L1=IGF1(I)                                                      
        L2=IGF2(I)                                                      
        L3=IGF3(I)                                                      
        ZXXX(I)=ZRB(L1,L2,L3)*DENOM                                     
   30 CONTINUE                                                          
C  -- Y -------------                                                   
          DO 1021 I=1,KSUM                                              
            ZV1D(I)=ZRY(I)*(RHV(I))/CHGB1(I)                        
 1021     CONTINUE                                                      
C*****---- FAST FOURIER TRANSFORMATION -----                            
      CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                  
     &          ,0,-1,1,IWORK,IERR)                                     
      IF (IERR.NE.0) THEN                                               
          WRITE (6,*) ' C3FFT(FFT)]  IERR = ',IERR                      
          STOP                                                          
      END IF                                                            
      DO 31 I=1,KG                                                      
        L1=IGF1(I)                                                      
        L2=IGF2(I)                                                      
        L3=IGF3(I)                                                      
        ZYYY(I)=ZRB(L1,L2,L3)*DENOM                                     
   31 CONTINUE                                                          
C  -- Z -------------                                                   
          DO 1022 I=1,KSUM                                              
            ZV1D(I)=ZRZ(I)*(RHV(I))/CHGB1(I)                        
 1022     CONTINUE                                                      
C*****---- FAST FOURIER TRANSFORMATION -----                            
      CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                  
     &          ,0,-1,1,IWORK,IERR)                                     
      IF (IERR.NE.0) THEN                                               
          WRITE (6,*) ' C3FFT(FFT)]  IERR = ',IERR                      
          STOP                                                          
      END IF                                                            
      DO 32 I=1,KG                                                      
        L1=IGF1(I)                                                      
        L2=IGF2(I)                                                      
        L3=IGF3(I)                                                      
        ZZZZ(I)=ZRB(L1,L2,L3)*DENOM                                     
   32 CONTINUE                                                          
      WRITE (6,*) 'ICOUNT = ',ICOUNT,' XSTPC END'                       
      RETURN                                                            
      END                                                               
C SUBROUTINE OCCUP2 ====*====3====*====4====*====5====*====6====*====7  
C                                                                       
C          1984.05.17 :   NORIAKI HAMADA                                
C                                                                       
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7  
C                                                                       
      SUBROUTINE FERMI(TOTCH,QWGT,WIDTH)
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
!XOCL SUBPROCESSOR PS(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IP=(PS,INDEX=1:KNV3,PART=BAND)
      DIMENSION EKK(KEG,KNV3),OCCUU(KEG,KNV3)
!XOCL LOCAL EKK(:,/IP),OCCUU(:,/IP)
      EQUIVALENCE (EKO,EKK),(OCCUP,OCCUU)
      DIMENSION QWGT(KNV3)                                              
      IF (ITER.EQ.1) WRITE(6,*) ' OCCUP2 :    WIDTH=',WIDTH             
      EMIN=EKO(NBD1,1)                                                  
      EMAX=EKO(NBD1,1)                                                  
!XOCL SPREAD DO /IP
      DO 10 K=1,KV3                                                     
      DO 10 I=NBD1,NBD2                                                 
        IF(EKK(I,K).LT.EMIN) EMIN=EKK(I,K)                              
        IF(EKK(I,K).GT.EMAX) EMAX=EKK(I,K)                              
   10 CONTINUE                                                          
!XOCL END SPREAD
      WSPIN  = 1.0D0                                                    
      EFERMI = EMAX                                                     
      E1     = EMIN                                                     
      E2     = EMAX                                                     
      JCOUNT = 1                                                        
C     WRITE(6,*) JCOUNT,TOT,EFERMI                                      
    1 CONTINUE                                                          
      TOT = 0.0D0                                                       
!XOCL SPREAD DO /IP
      DO 20 K=1,KV3                                                     
      DO 20 I=NBD1,NBD2                                                 
        E = EKK(I,K)                                                    
        CALL WIDTH2(E,EFERMI,WIDTH,DOS,WEIGHT)                          
        OCCUU(I,K) = (WEIGHT*WSPIN)*DBLE(KV3)*QWGT(K)                   
        TOT = TOT + 2.0D0*OCCUU(I,K)                                    
   20 CONTINUE                                                          
!XOCL END SPREAD SUM(TOT) 
      TOT=TOT/DFLOAT(KV3)                                               
C     WRITE(6,*) JCOUNT,TOT,EFERMI                                      
      IF(JCOUNT.EQ.1 .AND. TOT .LT. TOTCH) THEN                         
        WRITE(6,*) ' EMIN=',EMIN,'    EMAX=',EMAX                       
        WRITE(6,*) ' TOT=',TOT, '   TOTCH =', TOTCH                     
        STOP ' === STOP IN SUB.OCCUP2. (TOO FEW OF STATES) =='          
      END IF                                                            
C
      IF(DABS(TOT-TOTCH).LT.1.0D-10) GO TO 2                            
C                                     ----------------------->          
      IF(TOT.LT.TOTCH) THEN                                             
        E1     = EFERMI                                                 
        EFERMI = EFERMI + (E2-EFERMI)/2.0D0                             
      ELSE                                                              
        E2     = EFERMI                                                 
        EFERMI = EFERMI + (E1-EFERMI)/2.0D0                             
      END IF                                                            
      JCOUNT = JCOUNT + 1                                               
      GO TO 1                                                           
C                                     <----------------------           
    2 CONTINUE                                                          
C                                                                       
      EF = EFERMI                                                       
      WRITE( 6,100) EFERMI,TOT                                          
  100 FORMAT(1H ,'---------- THE FERMI ENERGY =',D23.16,F12.6)          
      RETURN                                                            
      END                                                               
C SUBROUTINE WIDTH2                                                     
C                                                                       
C          1983/05/18 :   NORIAKI HAMADA                                
C                                                                       
C----*----1----*----2----*----3----*----4----*----5----*----6----*----7 
C                                                                       
      SUBROUTINE WIDTH2                                                 
C INPUT                                                                 
     &                  ( E,EFF,W,                                      
C OUTPUT                                                                
     &                    DOS,OCC  )                                    
C                                                                       
C          EIGENENERGY IS BROADENED                                     
C                                                                       
C          E :  EIGENENERGY                                             
C          EF:  FERMI ENERGY                                            
C          W :  4*W IS BOTTOM LENGTH                                    
C          DOS    :  DENSITY OF STATES AT EF                            
C          OCCUP :  OCUPATION FRACTION OF ELECTRON                      
C                                                                       
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      IF(W .LE. 0.0D0) STOP ' === STOP AT SUB.WIDTHE. (W<=0.0D0) ==='   
      EE=(EFF-E)/W                                                      
      WW=4.0D0*W                                                        
      IF(EE .LE. -2.0D0) THEN                                           
        DOS    = 0.0D0                                                  
        OCC  = 0.0D0                                                    
      ELSE IF(EE .LT. -1.0D0) THEN                                      
        DOS    = (EE+2.0D0)**2/WW                                       
        OCC  = ((EE+2.0)**3)/(12.0D0)                                   
      ELSE IF(EE .LT. 1.0D0) THEN                                       
        DOS    = (2.0D0-EE**2)/WW                                       
        OCC  = (6.0D0+6.0D0*EE-(EE**3))/12.0D0                          
      ELSE IF(EE .LT. 2.0D0) THEN                                       
        DOS    = ((EE-2.0D0)**2)/WW                                     
        OCC  = (12.0D0+(EE-2.0D0)**3)/12.0D0                            
      ELSE                                                              
        DOS    = 0.0D0                                                  
        OCC  = 1.0D0                                                    
      END IF                                                            
      RETURN                                                            
      END                                                               
C     /////////////////////////////////////                             
C     // CALCULATION OF THE TOTAL ENERGY //                             
C     /////////////////////////////////////                             
      SUBROUTINE ENERGY(TOTCH,ETOT1,ETOT2,ETOTAL,SCHGPC,EPC)
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
!XOCL SUBPROCESSOR PS(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IP=(PS,INDEX=1:KNV3,PART=BAND)
      DIMENSION EKK(KEG,KNV3),OCCUU(KEG,KNV3)
!XOCL LOCAL EKK(:,/IP),OCCUU(:,/IP)
      EQUIVALENCE (EKO,EKK),(OCCUP,OCCUU)
C     /////////////////////////////////////////////////                 
C     // PSEUDOPOTENTIAL-CORRECTION AND EWALD ENERGY //                 
C     /////////////////////////////////////////////////                 
      ETOTAL = ETOT1*TOTCH + ETOT2                                      
C     ///////////////////////////////////////////////////////////       
C     // BAND ENERGY(1) <- "VXC" AND "VHR"   WITH "OLD CHARGE" //       
C     ///////////////////////////////////////////////////////////       
      CALL XCFFT(2,SCHGPC)                                 
      ETOTAL = ETOTAL - UNIVOL*DREAL(DCONJG(ZVXC(1))*ZCHG(1))           
     &          + UNIVOL*DREAL(DCONJG(ZEXC(1))*(ZCHG(1)+ZRHPC(1)))      
      DO 10 I=2,KG                                                      
      ETOTAL = ETOTAL - UNIVOL*                                         
     &         ( DREAL(DCONJG(ZVXC(I))*ZCHG(I)) +                       
     &           PAI4*REAL( (DCONJG(ZCHGO(I)))*ZCHG(I) )*RGG(I)  )    
     &                + UNIVOL*                                         
     &         ( DREAL(DCONJG(ZEXC(I))*(ZCHG(I)+ZRHPC(I))) +            
     &           PAI2*REAL( (DCONJG(ZCHG(I)))*ZCHG(I) )*RGG(I)  )     
10    CONTINUE                                                          
C     //////////////////////////////////////////////                    
C     // BAND ENERGY(3); BRILLOUIN ZONE SUMMATION //                    
C     //////////////////////////////////////////////                    
      FFF = FLOAT(KV3)
      TTT = 0.D0                                                  
!XOCL SPREAD DO /IP
      DO 2 I=1,KV3                                                
      DO 100 IBAN=NBD1,NBD2                                             
            TTT = TTT + OCCUU(IBAN,I)*                                  
     &      EKK(IBAN,I)                                                 
  100 CONTINUE                                                          
    2 CONTINUE                                                    
!XOCL END SPREAD SUM(TTT)
      ETOTAL = ETOTAL + 2.D0*TTT/FFF                                    
C FOR PARTIAL CORE CORRECTION     90.11.3                               
      ETOTAL = ETOTAL-EPC                                               
      WRITE(6,1104) ITER,ETOTAL,ETOTAL                                  
      WRITE(19,1104) ITER,ETOTAL,ETOTAL                                 
 1104 FORMAT(1H ,'TOTAL ENERGY FOR',I4,'-TH ITERATION=',                
     &            F12.7,5X,D15.7)                                       
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                  
      SUBROUTINE PSEUDO(TOTCH,ETOT1)                                    
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                  
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      DIMENSION VD(KIZA),VDNL(KIZA)                                     
C     DATA FOR ATOMS                                                    
      TOTCH = 0.D0                                                      
      ETOT1 = 0.D0                                                      
      UNIRV = 1.0D0/UNIVOL
      RG(1) =0.0D0                                                      
      RGG(1)=0.0D0                                                      
      DO 2000 I=2,KG                                                    
      RG(I) =1.0D0/(GR(I))                                       
      RGG(I)=1.0D0/(GR(I)*GR(I))                                 
 2000 CONTINUE                                                          
C                                                                       
      REWIND 16                                                         
      DO 500 II=1,KTYP                                                  
      ETOT1 = ETOT1 + FLOAT(IATOM(II))*PAI*ACHG(II)                     
     &      * ( AC(II,1)/BC(II,1) + AC(II,2)/BC(II,2) )/UNIVOL          
      TOTCH=TOTCH + FLOAT(IATOM(II))*ACHG(II)                           
C                                                                       
      IF (ITER.EQ.0) THEN                                               
        IF (NLSPD(II).EQ.1) THEN
          READ(16,*) MESHR,NMES,DX,RAD,VD,VDNL
          DO 1210 N=1,MESHR                                             
          VDD(N,II)=VD(N)                                               
          VDDNL(N,II)=VDNL(N)                                           
 1210 CONTINUE                                                          
        ELSE
          DO 1212 N=1,MESHR                                             
          VDD(N,II)  =0.0D0
          VDDNL(N,II)=0.0D0
 1212 CONTINUE                                                          
        END IF          
      ELSE                                                              
          DO 1211 N=1,MESHR                                             
          VD(N)=VDD(N,II)                                               
          VDNL(N)=VDDNL(N,II)                                           
 1211 CONTINUE                                                          
      END IF                                                            
C                                                                       
      S=0.0D0                                                           
      DO 1200 N=1,MESHR                                                 
      S=S + OMO(N)*VDNL(N)*(RAD(N)**3)                                  
 1200 CONTINUE                                                          
      PSC(1,II)=0.0D0                                                   
      DSC(1,II)=0.0D0                                                   
      ETOT1    =ETOT1+FLOAT(IATOM(II))*4.0D0*PAI*DX*S/UNIVOL            
C     CORE-PART NCP MATRIX ELEMENT; PSCO                                
          C1=AC(II,1)*4.D0*PAI*ACHG(II)                                 
          C2=AC(II,2)*4.D0*PAI*ACHG(II)                                 
          B1=-0.25D0/BC(II,1)                                           
          B2=-0.25D0/BC(II,2)                                           
      DO 82 I=2,KG                                                      
      GG=GR(I)**2                                                       
      PSC(I,II)=-( C1*EXP(B1*GG) + C2*EXP(B2*GG) )/(GG*UNIVOL)          
      DSC(I,II)=-( C1*(B1-RGG(I))*EXP(B1*GG)                     
     &            +C2*(B2-RGG(I))*EXP(B2*GG) )*RGG(I)*UNIRV            
82    CONTINUE                                                          
C                                                                       
      IF (NLSPD(II).EQ.1) THEN
C
      DO 1201 N=1,MESHR                                                 
      R=RAD(N)                                                          
      FAC1=OMO(N)*VDNL(N)*RAD(N)**3*4.0D0*PAI*DX/UNIVOL                 
      FAC2=OMO(N)*VDNL(N)*RAD(N)**4*2.0D0*PAI*DX                        
      DO 1110 I=1,KG                                                    
      XX(I)=R*GR(I)                                                      
 1110 CONTINUE                                                          
      CALL DNV(0,KG,XX,Y11)                                             
      CALL DNV(1,KG,XX,Y22)                                             
      DO 1120 I=2,KG                                                    
      PSC(I,II)=PSC(I,II)+FAC1*Y11(I)
      DSC(I,II)=DSC(I,II)-FAC2*Y22(I)*RG(I)*UNIRV                              
 1120 CONTINUE                                                          
 1201 CONTINUE                                                          
      END IF
  500 CONTINUE                                                          
C     ETOT1 -> TOTAL ENERGY CORRECTION FROM THE PSEUDOPOTENTIAL         
C     ETOT2 -> TOTAL ENERGY CONTRIBUTION FROM THE EWALD SUMMATION       
      WRITE(6,*) 'TOTAL VAL CHG PER UNIT-CELL    = ',TOTCH              
      WRITE(6,*) 'TOTAL ENERGY CORRECTION FROM PS= ',ETOT1              
      WRITE(6,*) '                         UNIVOL= ',UNIVOL             
      DO 75 I=1,KG                                                      
      ZPSCC(I) = DCMPLX(0.0D0,0.0D0)                                    
      ZDSCC(I) = DCMPLX(0.0D0,0.0D0)                                    
75    CONTINUE                                                          
      DO 73 IT=1,KTYP                                                   
      DO 72 I=1,KG                                                      
      ZPSCC(I) = ZPSCC(I) + PSC(I,IT        )*ZFM3(I,IT)                
      ZDSCC(I) = ZDSCC(I) + DSC(I,IT        )*ZFM3(I,IT)                
72    CONTINUE                                                          
73    CONTINUE                                                          
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                  
      SUBROUTINE PSELMD                                                 
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                  
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      DO 75 I=1,KG                                                      
      ZPSCC(I) = DCMPLX(0.0D0,0.0D0)                                    
      ZDSCC(I) = DCMPLX(0.0D0,0.0D0)
75    CONTINUE                                                          
      DO 73 IT=1,KTYP                                                   
      DO 72 I=1,KG                                                      
      ZPSCC(I) = ZPSCC(I) + PSC(I,IT        )*ZFM3(I,IT)                
      ZDSCC(I) = ZDSCC(I) + DSC(I,IT        )*ZFM3(I,IT)                
72    CONTINUE                                                          
73    CONTINUE                                                          
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                        
      SUBROUTINE DIAGON(IREC8)                                          
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                        
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      PARAMETER(KNMT=KNG11*(KNG11+1)/2,KNM9=9*KNG11)                    
!XOCL SUBPROCESSOR PS(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IP=(PS,INDEX=1:KNV3,PART=BAND)
      DIMENSION ZZ2(KNG1,KEG,KNV3),SSS(KNG1,KNV3,KTYP,10)
      DIMENSION EKK(KEG,KNV3),OCCUU(KEG,KNV3) 
!XOCL LOCAL ZZ2(:,:,/IP),SSS(:,/IP,:,:)
!XOCL LOCAL EKK(:,/IP),OCCUU(:,/IP)
      EQUIVALENCE (ZAJ,ZZ2),(SNL,SSS),(EKO,EKK)
     &           ,(OCCUP,OCCUU)
C     EIGEN-VALUE PROBLEM                                               
      LOGICAL*4 LW
      DIMENSION ZAAA(KNMT)                                              
     &         ,IFLG(KEG)
     &         ,ZVN(KNG11,KEG),ZZZ(KNG11,KNG11),EG(KEG)
     &         ,ZAWORK(KNG1)
     &         ,WWE(KNG11,7),LW(KNG11)
CEIG &         ,ZAWORK(KNG1)
CEIG &         ,WWE(KNG11,7),LW(KNG11)
CCHB &         ,WWW(3*KNG11),ZWW(5*KNG11),ZAWORK(KNG1)
CHEI &         ,EVR(KNG11,KEG),EVI(KNG11,KEG),EG(KEG)                   
CHEI &         ,WK1(KNG11,6),WK2(8*KEG),ZWK3(KNG11,2),IWK(6*KEG+KNG11)  
CHEI &         ,WOK(KNM9),CCC(KNG11,KNG11)                              
C                                                                      
      DIMENSION PPM(KNG11),IJG(KNG11)                                     
C                                                                       
      IF (KNG1.NE.KNG11) THEN                                           
          WRITE (6,*) 'WORNING]  KNG1 IS NOT EQUAL TO KNG11'            
      END IF                                                            
C                                                                       
      DO 4 I=1,KG                                                       
      ZCHGO(I) = ZCHG(I)                                                
      ZCHG(I)  = DCMPLX(0.0D0,0.0D0)                                    
    4 CONTINUE                                                          
C     ////////////////////////////                                      
C     // DIAGONALIZE AT K-POINT //                                      
C     ////////////////////////////
C     VPP-PARALLEL START
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 100 NNN=1,KNV3                                                 
C                                  IWRT(NNN) =NNN                      
                                   AKX = VX(NNN)                        
                                   AKY = VY(NNN)                        
                                   AKZ = VZ(NNN)                        
      IIKB = IBA(NNN)                                                   
      IIBA = IBA2(NNN)                                                  
C     PSEUDOPOTENTIAL NON-LOCAL PART MATRIX ELEMENT                     
      DO 71 I=1,KNG11*(KNG11+1)/2                                       
        ZAAA(I)=DCMPLX(0.0D0,0.0D0)                                     
 71   CONTINUE                                                          
      DO 272 I=1,KNG11                                                  
        DO 274 K=1,KEG                                                  
C         EVR(I,K)=0.0D0                                                
C         EVI(I,K)=0.0D0                                                
        ZVN(I,K)=CMPLX(0.0,0.0)
  274   CONTINUE                                                        
        DO 273 J=1,KNG11                                                
C         CCC(J,I)=0.0D0                                                
          ZZZ(J,I)=CMPLX(0.0,0.0)
  273   CONTINUE                                                        
  272 CONTINUE                                                          
C                                                                       
      DO 6 IT=1,KTYP                                                    
        CS=1.0D0/(WS(IT)*UNIVOL)                                        
        CP=1.0D0/(WP(IT)*UNIVOL)                                        
        CWL(1)=CS
        CWL(2)=CP
        CWL(3)=CP
        CWL(4)=CP
        IF (NLSPD(IT).EQ.2) THEN
            CD=1.0D0/(WD(IT)*UNIVOL)
            CWL(5)=CD
            CWL(6)=CD
            CWL(7)=CD
            CWL(8)=CD
            CWL(9)=CD
        END IF
          IF (NLSPD(IT).EQ.1) THEN
              LNUM = 4
          ELSE
              LNUM = 9
          END IF
        DO 5 I=1,IIBA                                                   
                       I1  = NBMAT(I,NNN)                               
                                           L1=IG1(I1)+KX1               
                                           L2=IG2(I1)+KY1               
                                           L3=IG3(I1)+KZ1               
          II=I*(I-1)/2                                                  
          DO 8700 L=1,LNUM
            DO 7600 J =1,I                                                
            PPMT   =   CWL(L)*SSS(I,NNN,IT,L)*SSS(J,NNN,IT,L)               
            ZAAA(II+J) =  ZAAA(II+J) + PPMT  *                          
     &                ZFM3(  IGPO( L1-IG1(NBMAT(J,NNN)),                
     &                             L2-IG2(NBMAT(J,NNN)),                
     &                             L3-IG3(NBMAT(J,NNN)) )  ,IT)         
 7600       CONTINUE                                                      
 8700     CONTINUE
    5   CONTINUE                                                        
    6 CONTINUE                                                          
C     MATRIX ELEMENT AND KINETIC ENERGY
      DO 230 I=1,IIBA                                                   
        PPM(I) = (   (AKX+GX(NBMAT(I,NNN)))**2                          
     &             + (AKY+GY(NBMAT(I,NNN)))**2                          
     &             + (AKZ+GZ(NBMAT(I,NNN)))**2  )/2.D0                  
  230 CONTINUE                                                          
      DO 222 I=1,IIBA                                                   
                           I1=NBMAT(I,NNN)                              
                           I2=I*(I-1)/2                                 
                                            L1=IG1(I1)+KX1              
                                            L2=IG2(I1)+KY1              
                                            L3=IG3(I1)+KZ1              
      ZAAA(I*(I+1)/2) = ZAAA(I*(I+1)/2) + PPM(I)                        
      DO 222 J=1,I                                                      
C     HARTREE, EXCHANGE AND CORE POTENTIAL                              
      IJG(J) = ( IGPO( L1-IG1(NBMAT(J,NNN)), L2-IG2(NBMAT(J,NNN)),      
     &                                       L3-IG3(NBMAT(J,NNN)) ) )   
      IF(IJG(J).NE.1) THEN                                              
      ZAAA(I2+J) = ZAAA(I2+J) + PAI4*ZCHGO(IJG(J))/GR(IJG(J))**2        
      END IF                                                            
      ZAAA(I2+J) = ZAAA(I2+J) + ZVXC(IJG(J)) + ZPSCC(IJG(J))            
222   CONTINUE                                                          
C=====DIAGONALIZATION SSL2                                              
      DO 7701 I=1,IIBA                                                  
      MM = I*(I-1)/2                                                    
      DO 7702 J=1,I                                                     
      ZZZ(I,J)=ZAAA(MM+J)
      IF (I.NE.J) ZZZ(J,I)=CONJG(ZZZ(I,J))
C     RMAT    = DREAL(ZAAA(MM+J))                                       
C     CMAT    = DIMAG(ZAAA(MM+J))                                       
C     CCC(I,J)  = RMAT                                                  
C     IF (I.NE.J) CCC(J,I) = CMAT                                       
 7702 CONTINUE                                                          
 7701 CONTINUE                                                          
C""""""""""""""""""""""""""""""""""""""""""""""""""""""""""             
      KBAS = KNG11                                                      
      KGB  = KEG                                                        
      IOPT = 2                                                          
      EPS  = 1.0D-15                                                    
C     CALL HZES1M(ZAAA,IIBA,KGB,KGB,EPS,IOPT,EG,EVR,KBAS,EVI,IFLG       
C    &           ,WK1,WK2,ZWK3,IWK,ICON)                                
C     WRITE(6,*) KBAS,IIBA,KGB                                          
C     CALL DHEIG2(CCC,KBAS,IIBA,-KGB,EG,EVR,EVI,WOK,ICON)               
C     CALL CHOBSD(ZZZ,KBAS,IIBA,EG,-KGB,ZVN,KGB,EPS,WWW,ZWW,ICON)
      CALL EIGCH(ZZZ,IIBA,KNG11,-KGB,KGB,EPS,WWE,LW,EG,ZVN,ICON)
C""""""""""""""""""""""""""""""""""""""""""""""""""""""""""             
      IF(ICON.NE.0) THEN                                                
      WRITE(6,*) '*****!! AT HZES1M,  ICON=',ICON                       
      END IF                                                            
C                                                                       
      DO 7164 M=1,KEG                                                   
      EKK(M,NNN)=EG(M)                                                  
C      WRITE(6,*) 'NNN = ',NNN,'IBAN = ',M                               
C      WRITE(6,*) 'ENERGY = ',EKK(M,NNN)
 7164 CONTINUE                                                          
C     OUTPUT EIGENVECTORS ON FILE 80 (DIRECT-ACCESS) (KNG11 <---> IIBA) 
C     WRITE(6,*) 'ZAJ = '                                               
      DO 231 IBAN = NBD1,NBD2                                           
        DO 33 J=1,IIBA                                                  
C         ZAJ(J,IBAN,NNN) = EVR(J,IBAN) + ZI*EVI(J,IBAN)                
          ZZ2(J,IBAN,NNN) = ZVN(J,IBAN)
C*********ZV1(J)          = EVR(J,IBAN) + ZI*EVI(J,IBAN)                
   33   CONTINUE                                                        
        DO 97 M=1,IIKB                                                  
          ZAWORK(M)=0.D0                                                
   97   CONTINUE                                                        
        DO 91 J=1,IIBA                                                  
          DO 92 M=1,IIKB                                                
            IF (NBASE(M,NNN) .EQ. NBMAT(J,NNN)) THEN                    
              ZAWORK(M)=ZZ2(J,IBAN,NNN)                                 
            ENDIF                                                       
 92       CONTINUE                                                      
 91     CONTINUE                                                        
        DO 93 M=1,IIKB                                                  
          ZZ2(M,IBAN,NNN)=ZAWORK(M)                                     
 93     CONTINUE                                                        
C**************DO 34 J=IIBA+1,KNG1                                      
C                ZV1(J)          = DCMPLX(0.0D0,0.0D0)                  
C  34          CONTINUE                                                 
C              IREC = KV3*(IBAN-1)+IWRT(NNN)                            
C              IF (IREC8.NE.0) IREC=IREC8*IREC-IREC8+1                  
C**************WRITE(80,REC=IREC)  ZV1                                  
  231 CONTINUE                                                          
  100 CONTINUE                                                          
!XOCL END SPREAD
C!XOCL END PARALLEL
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      SUBROUTINE KBMAT                                                  
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                                    
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
C     VPP-PARALLEL
!XOCL SUBPROCESSOR PS(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IP=(PS,INDEX=1:KNV3,PART=BAND)
      DIMENSION SSS(KNG1,KNV3,KTYP,10),RAA(KNG1,KNV3)
!XOCL LOCAL SSS(:,/IP,:,:),RAA(:,/IP)
      EQUIVALENCE (SNL,SSS),(RAK,RAA)
C     K-B TYPE ADDITION                                                 
      DIMENSION RKS(KIZA),RKP(KIZA),RKD(KIZA),DVS(KIZA)
     &         ,DVP(KIZA),DVD(KIZA)
      NSEK=2
      SC=DSQRT(PAI4)                                                    
      PC=DSQRT(3.0D0*PAI4)                                              
      DC=DSQRT(5.0D0*PAI4)
      CO1=DSQRT(3.0D0)
      REWIND 15                                                         
C     INTEGRATION OF K-B TYPE POTENTIAL                                 
      DO 100 ITY = 1,KTYP                                               
        IF (NLSPD(ITY).EQ.1) THEN
          READ (15,*) MESHR,NMES,DX,WS(ITY),WP(ITY)
     &             ,RAD,DVS,DVP
          WD(ITY)=0.0D0
          DO 111 N=1,MESHR                                                
            WVS(N,ITY)=DVS(N)                                             
            WVP(N,ITY)=DVP(N)                                             
            WVD(N,ITY)=0.0D0
  111 CONTINUE
        ELSE
          READ (15,*) MESHR,NMES,DX,WS(ITY),WP(ITY),WD(ITY)
     &               ,RAD,DVS,DVP,DVD
          DO 110 N=1,MESHR                                                
            WVS(N,ITY)=DVS(N)                                             
            WVP(N,ITY)=DVP(N)                                             
            WVD(N,ITY)=DVD(N)
  110 CONTINUE                                                        
        END IF
  100 CONTINUE                                                          
C
!XOCL SPREAD DO /IP      
      DO 1003 IK=1,KV3                                                  
      AKX = VX(IK)                                                      
      AKY = VY(IK)                                                      
      AKZ = VZ(IK)                                                      
C      IIBA = IBA(IK)
      DO 1005 I=1,IBA2(IK)
      RK=SQRT((AKX+GX(NBMAT(I,IK)))**2                                  
     &   +(AKY+GY(NBMAT(I,IK)))**2 + (AKZ+GZ(NBMAT(I,IK)))**2 )     
*VOCL STMT,IF(90)
      IF (RK.NE.0.0D0) THEN                                             
          RAA(I,IK)=1.0D0/RK
      ELSE                                                              
          RAA(I,IK)=0.0D0                                               
      END IF                                                            
 1005 CONTINUE                                                          
 1003 CONTINUE                                                          
!XOCL END SPREAD
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 200 NNN=1,KV3                                                  
                                   AKX  = VX(NNN)                       
                                   AKY  = VY(NNN)                       
                                   AKZ  = VZ(NNN)                       
                                   IIBA = IBA2(NNN)                     
        DO 1000 J =1,IIBA                                               
          QX(J)  = AKX + GX(NBMAT(J,NNN))                               
          QY(J)  = AKY + GY(NBMAT(J,NNN))                               
          QZ(J)  = AKZ + GZ(NBMAT(J,NNN))                               
          QA1(J) = 0.5D0*(2.0D0*QZ(J)*QZ(J)-QX(J)*QX(J)-QY(J)*QY(J))
          QA2(J) = 0.5D0*CO1*(QX(J)*QX(J)-QY(J)*QY(J))
          QA3(J) = CO1*QX(J)*QY(J)
          QA4(J) = CO1*QX(J)*QZ(J)
          QA5(J) = CO1*QY(J)*QZ(J)
          AK2(J) = SQRT( QX(J)**2 + QY(J)**2 + QZ(J)**2 )               
 1000   CONTINUE                                                        
        DO 1100 ITY=1,KTYP                                              
          DO 1111 KN=1,KNG11                                            
            SSS(KN,NNN,ITY,1)=0.0D0                                     
            SSS(KN,NNN,ITY,2)=0.0D0                                     
            SSS(KN,NNN,ITY,3)=0.0D0                                     
            SSS(KN,NNN,ITY,4)=0.0D0                                     
 1111     CONTINUE                                                      
C                                                                       
      DO 7320 N=1,MESHR                                                 
      R=RAD(N)                                                          
      FACS= OMO(N)*WVS(N,ITY)*R**NSEK*SC*DX                             
      FACP= OMO(N)*WVP(N,ITY)*R**NSEK*PC*DX                             
      DO 7777 KN=1,IIBA                                                 
        X(KN)=AK2(KN)*R                                                 
 7777 CONTINUE                                                          
      CALL DNV(0,IIBA,X,Y1)                                           
      CALL DNV(1,IIBA,X,Y2)                                           
      DO 7330 KN=1,IIBA                                                 
C     RKS(N)=$DMSJM(0,AK2(KN)*RAD(N))                                   
      SSS(KN,NNN,ITY,1) =SSS(KN,NNN,ITY,1) +FACS*Y1(KN)                 
      SSS(KN,NNN,ITY,2) =SSS(KN,NNN,ITY,2)                              
     &                  +FACP*Y2(KN)*QX(KN)*RAA(KN,NNN)                 
      SSS(KN,NNN,ITY,3) =SSS(KN,NNN,ITY,3)                              
     &                  +FACP*Y2(KN)*QY(KN)*RAA(KN,NNN)                 
      SSS(KN,NNN,ITY,4) =SSS(KN,NNN,ITY,4)                              
     &                  +FACP*Y2(KN)*QZ(KN)*RAA(KN,NNN)                 
 7330 CONTINUE                                                          
 7320 CONTINUE                                                          
C
      IF (NLSPD(ITY).EQ.2) THEN
          DO 1112 KN=1,KNG11                                            
            SSS(KN,NNN,ITY,5)=0.0D0                                     
            SSS(KN,NNN,ITY,6)=0.0D0                                     
            SSS(KN,NNN,ITY,7)=0.0D0                                     
            SSS(KN,NNN,ITY,8)=0.0D0                                     
            SSS(KN,NNN,ITY,9)=0.0D0                                     
 1112    CONTINUE                                                      
C                                                                       
      DO 7321 N=1,MESHR                                                 
      R=RAD(N)                                                          
      FACD= OMO(N)*WVD(N,ITY)*R**NSEK*DC*DX
      DO 7778 KN=1,IIBA                                                 
        X(KN)=AK2(KN)*R                                                 
 7778 CONTINUE                                                          
      CALL DNV(2,IIBA,X,Y3)                                           
      DO 7331 KN=1,IIBA                                                 
C     RKS(N)=$DMSJM(0,AK2(KN)*RAD(N))                                   
      SSS(KN,NNN,ITY,5) =SSS(KN,NNN,ITY,5)
     &                  +FACD*Y3(KN)*QA1(KN)*RAA(KN,NNN)*RAA(KN,NNN)
      SSS(KN,NNN,ITY,6) =SSS(KN,NNN,ITY,6)
     &                  +FACD*Y3(KN)*QA2(KN)*RAA(KN,NNN)*RAA(KN,NNN)
      SSS(KN,NNN,ITY,7) =SSS(KN,NNN,ITY,7)
     &                  +FACD*Y3(KN)*QA3(KN)*RAA(KN,NNN)*RAA(KN,NNN)
      SSS(KN,NNN,ITY,8) =SSS(KN,NNN,ITY,8)
     &                  +FACD*Y3(KN)*QA4(KN)*RAA(KN,NNN)*RAA(KN,NNN)
      SSS(KN,NNN,ITY,9) =SSS(KN,NNN,ITY,9)
     &                  +FACD*Y3(KN)*QA5(KN)*RAA(KN,NNN)*RAA(KN,NNN)    
 7331 CONTINUE                                                          
 7321 CONTINUE                                                          
      END IF
 1100 CONTINUE
  200 CONTINUE                                                          
!XOCL END SPREAD
C!XOCL END PARALLEL
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      SUBROUTINE KBINT
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
C     K-B TYPE ADDITION                                                 
C     VPP-PARALLEL
!XOCL SUBPROCESSOR PS(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IP=(PS,INDEX=1:KNV3,PART=BAND)
      DIMENSION SSS(KNG1,KNV3,KTYP,10),SS2(KNG1,KNV3,KTYP,9)
     &         ,SS3(KNG1,KNV3,KTYP,4)
     &         ,RAA(KNG1,KNV3)
!XOCL LOCAL SSS(:,/IP,:,:),SS2(:,/IP,:,:),SS3(:,/IP,:)
!XOCL LOCAL RAA(:,/IP)
      EQUIVALENCE (SNL,SSS),(SNL2,SS2),(SNL3,SS3)
     &           ,(RAK,RAA)
C
C     K-B TYPE ADDITION                                                 
C     INTEGRATION OF K-B TYPE POTENTIAL                                 
C     MODIFIED TO STRESS CALCULATION FOR KB-SEPARABLE FORM              
      NSEK=2                                                            
      CO0=DSQRT(2.0D0) 
      CO1=DSQRT(3.0D0)
      SC=DSQRT(PAI4)                                                    
      PC=DSQRT(3.0D0*PAI4)                                              
      DC=DSQRT(5.0D0*PAI4)
C
!XOCL SPREAD DO /IP      
      DO 1003 IK=1,KV3                                                  
      AKX = VX(IK)                                                      
      AKY = VY(IK)                                                      
      AKZ = VZ(IK)                                                      
C      IIBA = IBA(IK)
      DO 1005 I=1,IBA(IK)
      RK=SQRT((AKX+GX(NBASE(I,IK)))**2                                  
     &       +(AKY+GY(NBASE(I,IK)))**2 + (AKZ+GZ(NBASE(I,IK)))**2 )     
*VOCL STMT,IF(90)
      IF (RK.NE.0.0D0) THEN                                             
          RAA(I,IK)=1.0D0/RK
      ELSE                                                              
          RAA(I,IK)=0.0D0                                               
      END IF                                                            
 1005 CONTINUE                                                          
 1003 CONTINUE                                                          
!XOCL END SPREAD
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 1006 NNN=1,KV3                                                 
      AKX = VX(NNN)                                                     
      AKY = VY(NNN)                                                     
      AKZ = VZ(NNN)                                                     
      IIBA = IBA(NNN)                                                   
C
        DO 7500 J =1,IIBA                                               
        QX(J)=AKX+GX(NBASE(J,NNN))                                      
        QY(J)=AKY+GY(NBASE(J,NNN))                                      
        QZ(J)=AKZ+GZ(NBASE(J,NNN))                                      
      QA1(J)=0.5D0*(2.0D0*QZ(J)*QZ(J)-QX(J)*QX(J)-QY(J)*QY(J))
      QA2(J)=0.5D0*CO1*(QX(J)*QX(J)-QY(J)*QY(J))
      QA3(J)=CO1*QX(J)*QY(J)
      QA4(J)=CO1*QX(J)*QZ(J)
      QA5(J)=CO1*QY(J)*QZ(J)
        AK2(J) = SQRT( QX(J)**2 + QY(J)**2 + QZ(J)**2 )               
 7500 CONTINUE                                                          
C      REWIND 15                                                        
      DO 7300 ITY=1,KTYP                                                
C     ATTN]   IL=1,4 (S,P ONLY)                                         
      DO 7311 KN=1,KNG1                                                 
      SSS(KN,NNN,ITY,1)=0.0D0                                           
      SSS(KN,NNN,ITY,2)=0.0D0                                           
      SSS(KN,NNN,ITY,3)=0.0D0                                           
      SSS(KN,NNN,ITY,4)=0.0D0                                           
      SS2(KN,NNN,ITY,1)=0.0D0                                          
      SS2(KN,NNN,ITY,2)=0.0D0                                          
      SS2(KN,NNN,ITY,3)=0.0D0                                          
      SS2(KN,NNN,ITY,4)=0.0D0                                          
      SS3(KN,NNN,ITY,1)=0.0D0                                          
 7311 CONTINUE                                                          
C                                                                       
      DO 7320 N=1,MESHR                                                 
      R=RAD(N)                                                          
      FACS= OMO(N)*WVS(N,ITY)*R**NSEK*SC*DX                             
      FACP= OMO(N)*WVP(N,ITY)*R**NSEK*PC*DX                             
      FAC1=-OMO(N)*WVS(N,ITY)*R**(NSEK+1)*SC*DX                         
      FAC2= OMO(N)*WVP(N,ITY)*R**(NSEK+1)*PC*DX                         
      DO 7777 KN=1,IIBA                                                 
        X(KN)=AK2(KN)*R                                                 
 7777 CONTINUE                                                          
      CALL DNV(0,IIBA,X,Y1)                                           
      CALL DNV(1,IIBA,X,Y2)                                           
      CALL DNV(2,IIBA,X,Y3)                                           
      DO 7878 KN=1,IIBA                                                 
      Y3(KN)=(Y1(KN)-2.0D0*Y3(KN))/3.0D0                                
 7878 CONTINUE                                                          
      DO 7330 KN=1,IIBA                                                 
C     RKS(N)=$DMSJM(0,AK2(KN)*RAD(N))                                   
      SSS(KN,NNN,ITY,1) =SSS(KN,NNN,ITY,1) +FACS*Y1(KN)                 
      SSS(KN,NNN,ITY,2) =SSS(KN,NNN,ITY,2)                              
     &                  +FACP*Y2(KN)*QX(KN)*RAA(KN,NNN)                 
      SSS(KN,NNN,ITY,3) =SSS(KN,NNN,ITY,3)                              
     &                  +FACP*Y2(KN)*QY(KN)*RAA(KN,NNN)                 
      SSS(KN,NNN,ITY,4) =SSS(KN,NNN,ITY,4)                              
     &                  +FACP*Y2(KN)*QZ(KN)*RAA(KN,NNN)                 
      SS2(KN,NNN,ITY,1) =SS2(KN,NNN,ITY,1)+FAC1*Y2(KN)                 
      SS2(KN,NNN,ITY,2) =SS2(KN,NNN,ITY,2)                             
     &                  +FAC2*Y3(KN)*QX(KN)*RAA(KN,NNN)                 
      SS2(KN,NNN,ITY,3) =SS2(KN,NNN,ITY,3)                             
     &                  +FAC2*Y3(KN)*QY(KN)*RAA(KN,NNN)                 
      SS2(KN,NNN,ITY,4) =SS2(KN,NNN,ITY,4)                             
     &                  +FAC2*Y3(KN)*QZ(KN)*RAA(KN,NNN)                 
      SS3(KN,NNN,ITY,1) =SS3(KN,NNN,ITY,1)+FACP*Y2(KN)                 
 7330 CONTINUE                                                          
 7320 CONTINUE                                                          
C
      IF (NLSPD(ITY).EQ.2) THEN
C
      DO 7314 KN=1,KNG1                                                 
      SSS(KN,NNN,ITY,5)=0.0D0
      SSS(KN,NNN,ITY,6)=0.0D0
      SSS(KN,NNN,ITY,7)=0.0D0
      SSS(KN,NNN,ITY,8)=0.0D0
      SSS(KN,NNN,ITY,9)=0.0D0
      SSS(KN,NNN,ITY,10)=0.0D0
      SS2(KN,NNN,ITY,5)=0.0D0
      SS2(KN,NNN,ITY,6)=0.0D0
      SS2(KN,NNN,ITY,7)=0.0D0
      SS2(KN,NNN,ITY,8)=0.0D0
      SS2(KN,NNN,ITY,9)=0.0D0
      SS3(KN,NNN,ITY,2)=0.0D0
      SS3(KN,NNN,ITY,3)=0.0D0
      SS3(KN,NNN,ITY,4)=0.0D0
 7314 CONTINUE                                                          
C                                                                       
      DO 7321 N=1,MESHR                                                 
      R=RAD(N)                                                          
      FACS= OMO(N)*WVS(N,ITY)*R**NSEK*SC*DX                             
      FACP= OMO(N)*WVP(N,ITY)*R**NSEK*PC*DX                             
      FACD= OMO(N)*WVD(N,ITY)*R**NSEK*DC*DX
      FAC1=-OMO(N)*WVS(N,ITY)*R**(NSEK+1)*SC*DX                         
      FAC2= OMO(N)*WVP(N,ITY)*R**(NSEK+1)*PC*DX                         
      FAC3= OMO(N)*WVD(N,ITY)*R**(NSEK+1)*DC*DX 
      DO 8778 KN=1,IIBA                                                 
        X(KN)=AK2(KN)*R                                                 
 8778 CONTINUE                                                          
      CALL DNV(0,IIBA,X,Y1)                                           
      CALL DNV(1,IIBA,X,Y2)                                           
      CALL DNV(2,IIBA,X,Y3)                                           
      CALL DNV(3,IIBA,X,Y4)
      DO 7788 KN=1,IIBA                                                 
      YD(KN)= Y3(KN)
      Y3(KN)=(Y1(KN)-2.0D0*Y3(KN))/3.0D0                                
      Y4(KN)=(2.0D0*Y2(KN)-3.0D0*Y4(KN))/5.0D0
 7788 CONTINUE                                                          
      DO 7331 KN=1,IIBA                                                 
C     ATTN] FOR SSL2 (IN ISSP)                                          
      SSS(KN,NNN,ITY,5) =SSS(KN,NNN,ITY,5)
     &                  +FACD*YD(KN)*QA1(KN)*RAA(KN,NNN)*RAA(KN,NNN)
      SSS(KN,NNN,ITY,6) =SSS(KN,NNN,ITY,6)
     &                  +FACD*YD(KN)*QA2(KN)*RAA(KN,NNN)*RAA(KN,NNN)
      SSS(KN,NNN,ITY,7) =SSS(KN,NNN,ITY,7)
     &                  +FACD*YD(KN)*QA3(KN)*RAA(KN,NNN)*RAA(KN,NNN)
      SSS(KN,NNN,ITY,8) =SSS(KN,NNN,ITY,8)
     &                  +FACD*YD(KN)*QA4(KN)*RAA(KN,NNN)*RAA(KN,NNN)
      SSS(KN,NNN,ITY,9) =SSS(KN,NNN,ITY,9)
     &                  +FACD*YD(KN)*QA5(KN)*RAA(KN,NNN)*RAA(KN,NNN)    
      SSS(KN,NNN,ITY,10)=SSS(KN,NNN,ITY,10)+FACD*YD(KN)  
      SS2(KN,NNN,ITY,5)=SS2(KN,NNN,ITY,5)
     &                  +FAC3*Y4(KN)*QA1(KN)*RAA(KN,NNN)*RAA(KN,NNN)
      SS2(KN,NNN,ITY,6)=SS2(KN,NNN,ITY,6)
     &                  +FAC3*Y4(KN)*QA2(KN)*RAA(KN,NNN)*RAA(KN,NNN)
      SS2(KN,NNN,ITY,7)=SS2(KN,NNN,ITY,7)
     &                  +FAC3*Y4(KN)*QA3(KN)*RAA(KN,NNN)*RAA(KN,NNN)
      SS2(KN,NNN,ITY,8)=SS2(KN,NNN,ITY,8)
     &                  +FAC3*Y4(KN)*QA4(KN)*RAA(KN,NNN)*RAA(KN,NNN)
      SS2(KN,NNN,ITY,9)=SS2(KN,NNN,ITY,9)
     &                  +FAC3*Y4(KN)*QA5(KN)*RAA(KN,NNN)*RAA(KN,NNN)
      SS3(KN,NNN,ITY,2)=SS3(KN,NNN,ITY,2)+CO1*FACD*QX(KN)*YD(KN)
     &                  *RAA(KN,NNN)
      SS3(KN,NNN,ITY,3)=SS3(KN,NNN,ITY,3)+CO1*FACD*QY(KN)*YD(KN)
     &                  *RAA(KN,NNN)
      SS3(KN,NNN,ITY,4)=SS3(KN,NNN,ITY,4)+CO1*FACD*QZ(KN)*YD(KN)
     &                  *RAA(KN,NNN)
 7331 CONTINUE                                                          
 7321 CONTINUE                                                          
      END IF
 7300 CONTINUE                                                          
 1006 CONTINUE
!XOCL END SPREAD
C!XOCL END PARALLEL
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                                   
      SUBROUTINE FORCE(IREC8)                                           
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                                   
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
C     VPP-PARALLEL
!XOCL SUBPROCESSOR PS(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IP=(PS,INDEX=1:KNV3,PART=BAND)
      DIMENSION SSS(KNG1,KNV3,KTYP,10)
     &         ,ZZZ(KNG1,KEG,KNV3),ZFC(KEG,KNV3,KATM,10)
     &         ,OCCUU(KEG,KNV3)  
!XOCL LOCAL SSS(:,/IP,:,:),ZZZ(:,:,/IP),ZFC(:,/IP,:,:)
!XOCL LOCAL OCCUU(:,/IP)
      EQUIVALENCE (SNL,SSS),(ZAJ,ZZZ),(ZFBB,ZFC)
     &           ,(OCCUP,OCCUU)
      DIMENSION ZFX1(KEG),ZFX2(KEG),ZFX3(KEG),ZFX4(KEG)                 
      DIMENSION ZFY1(KEG),ZFY2(KEG),ZFY3(KEG),ZFY4(KEG)                 
      DIMENSION ZFZ1(KEG),ZFZ2(KEG),ZFZ3(KEG),ZFZ4(KEG)                 
      DIMENSION ZFB1(KEG),ZFB2(KEG),ZFB3(KEG),ZFB4(KEG)                 
C
      DIMENSION ZFX5(KEG),ZFX6(KEG),ZFX7(KEG),ZFX8(KEG),ZFX9(KEG)
      DIMENSION ZFY5(KEG),ZFY6(KEG),ZFY7(KEG),ZFY8(KEG),ZFY9(KEG)
      DIMENSION ZFZ5(KEG),ZFZ6(KEG),ZFZ7(KEG),ZFZ8(KEG),ZFZ9(KEG)
      DIMENSION ZFB5(KEG),ZFB6(KEG),ZFB7(KEG),ZFB8(KEG),ZFB9(KEG)
      DIMENSION ZFBA(KEG)
C                                                                       
      DO 1000 IA=1,KATM                                                 
        ZFORC2(IA,1)=DCMPLX(0.0D0,0.0D0)                                
        ZFORC2(IA,2)=DCMPLX(0.0D0,0.0D0)                                
        ZFORC2(IA,3)=DCMPLX(0.0D0,0.0D0)                                
 1000 CONTINUE                                                          
C     VPP-PAEALLEL START
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 2100 IK=1,KV3                                                  
        IIBA =  IBA(IK)                                                 
C*******DO 1 IBAN=NBD1,NBD2                                            
C         IREC = KV3*(IBAN-1)+IK                                        
C         IF (IREC8.NE.0) IREC=IREC8*IREC-IREC8+1                       
C         READ(80,REC=IREC) ZV1                                         
C         DO 2 I=1,IIBA                                                 
C           ZEVC(I,IBAN) = ZV1(I)                                       
C   2     CONTINUE                                                      
C***1   CONTINUE                                                        
        DO 2110 IA=1,KATM                                               
          CS=1.0D0/(WS(KFTYPE(IA))*UNIVOL)                              
          CP=1.0D0/(WP(KFTYPE(IA))*UNIVOL)                              
          DO 2115 IBAN=NBD1,NBD2                                        
            ZFX1(IBAN)=DCMPLX(0.0D0,0.0D0)                              
            ZFX2(IBAN)=DCMPLX(0.0D0,0.0D0)                              
            ZFX3(IBAN)=DCMPLX(0.0D0,0.0D0)                              
            ZFX4(IBAN)=DCMPLX(0.0D0,0.0D0)                              
            ZFY1(IBAN)=DCMPLX(0.0D0,0.0D0)                             
            ZFY2(IBAN)=DCMPLX(0.0D0,0.0D0)                              
            ZFY3(IBAN)=DCMPLX(0.0D0,0.0D0)                              
            ZFY4(IBAN)=DCMPLX(0.0D0,0.0D0)                              
            ZFZ1(IBAN)=DCMPLX(0.0D0,0.0D0)                              
            ZFZ2(IBAN)=DCMPLX(0.0D0,0.0D0)                              
            ZFZ3(IBAN)=DCMPLX(0.0D0,0.0D0)                              
            ZFZ4(IBAN)=DCMPLX(0.0D0,0.0D0)                              
            ZFB1(IBAN)=DCMPLX(0.0D0,0.0D0)                              
            ZFB2(IBAN)=DCMPLX(0.0D0,0.0D0)                              
            ZFB3(IBAN)=DCMPLX(0.0D0,0.0D0)                              
            ZFB4(IBAN)=DCMPLX(0.0D0,0.0D0)                              
 2115     CONTINUE                                                      
          DO 2200 IBAN=NBD1,NBD2                                        
            DO 1510 I=1,IIBA                                            
              I1=NBASE(I,IK)                                            
              ZTMP=ZZZ(I,IBAN,IK)*DCONJG( ZFM2( I1 ,IA ) )             
              ZTMP1=ZTMP*SSS(I,IK,KFTYPE(IA),1)                         
              ZTMP2=ZTMP*SSS(I,IK,KFTYPE(IA),2)                         
              ZTMP3=ZTMP*SSS(I,IK,KFTYPE(IA),3)                         
              ZTMP4=ZTMP*SSS(I,IK,KFTYPE(IA),4)                         
              ZFB1(IBAN)=ZFB1(IBAN)+ZTMP1                               
              ZFB2(IBAN)=ZFB2(IBAN)+ZTMP2                               
              ZFB3(IBAN)=ZFB3(IBAN)+ZTMP3                               
              ZFB4(IBAN)=ZFB4(IBAN)+ZTMP4                               
              ZFX1(IBAN)=ZFX1(IBAN)+GX(I1)*ZTMP1                        
              ZFX2(IBAN)=ZFX2(IBAN)+GX(I1)*ZTMP2                        
              ZFX3(IBAN)=ZFX3(IBAN)+GX(I1)*ZTMP3                        
              ZFX4(IBAN)=ZFX4(IBAN)+GX(I1)*ZTMP4                        
              ZFY1(IBAN)=ZFY1(IBAN)+GY(I1)*ZTMP1                        
              ZFY2(IBAN)=ZFY2(IBAN)+GY(I1)*ZTMP2                        
              ZFY3(IBAN)=ZFY3(IBAN)+GY(I1)*ZTMP3                        
              ZFY4(IBAN)=ZFY4(IBAN)+GY(I1)*ZTMP4                        
              ZFZ1(IBAN)=ZFZ1(IBAN)+GZ(I1)*ZTMP1                        
              ZFZ2(IBAN)=ZFZ2(IBAN)+GZ(I1)*ZTMP2                        
              ZFZ3(IBAN)=ZFZ3(IBAN)+GZ(I1)*ZTMP3                        
              ZFZ4(IBAN)=ZFZ4(IBAN)+GZ(I1)*ZTMP4                        
 1510       CONTINUE                                                    
 2200     CONTINUE                                                      
          DO 2205 IBAN=NBD1,NBD2                                        
            ZFC(IBAN,IK,IA,1)=CS*ZFB1(IBAN)                             
            ZFC(IBAN,IK,IA,2)=CP*ZFB2(IBAN)                             
            ZFC(IBAN,IK,IA,3)=CP*ZFB3(IBAN)                             
            ZFC(IBAN,IK,IA,4)=CP*ZFB4(IBAN)                             
            ZFORC2(IA,1)=ZFORC2(IA,1)+OCCUU(IBAN,IK)*                   
     &                   ( DCONJG(ZFX1(IBAN))*ZFC(IBAN,IK,IA,1)         
     &                    +DCONJG(ZFX2(IBAN))*ZFC(IBAN,IK,IA,2)         
     &                    +DCONJG(ZFX3(IBAN))*ZFC(IBAN,IK,IA,3)         
     &                    +DCONJG(ZFX4(IBAN))*ZFC(IBAN,IK,IA,4)  )      
            ZFORC2(IA,2)=ZFORC2(IA,2)+OCCUU(IBAN,IK)*                   
     &                   ( DCONJG(ZFY1(IBAN))*ZFC(IBAN,IK,IA,1)         
     &                    +DCONJG(ZFY2(IBAN))*ZFC(IBAN,IK,IA,2)         
     &                    +DCONJG(ZFY3(IBAN))*ZFC(IBAN,IK,IA,3)         
     &                    +DCONJG(ZFY4(IBAN))*ZFC(IBAN,IK,IA,4)  )      
            ZFORC2(IA,3)=ZFORC2(IA,3)+OCCUU(IBAN,IK)*                   
     &                   ( DCONJG(ZFZ1(IBAN))*ZFC(IBAN,IK,IA,1)         
     &                    +DCONJG(ZFZ2(IBAN))*ZFC(IBAN,IK,IA,2)         
     &                    +DCONJG(ZFZ3(IBAN))*ZFC(IBAN,IK,IA,3)         
     &                    +DCONJG(ZFZ4(IBAN))*ZFC(IBAN,IK,IA,4)  )      
 2205     CONTINUE                                                      
C
        IF (NLSPD(KFTYPE(IA)).EQ.2) THEN
C
          CD=1.0D0/(WD(KFTYPE(IA))*UNIVOL)
          DO 2135 IBAN=NBD1,NBD2                                        
            ZFX5(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFX6(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFX7(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFX8(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFX9(IBAN)=DCMPLX(0.0D0,0.0D0) 
            ZFY5(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFY6(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFY7(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFY8(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFY9(IBAN)=DCMPLX(0.0D0,0.0D0) 
            ZFZ5(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFZ6(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFZ7(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFZ8(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFZ9(IBAN)=DCMPLX(0.0D0,0.0D0) 
            ZFB5(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFB6(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFB7(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFB8(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFB9(IBAN)=DCMPLX(0.0D0,0.0D0)
            ZFBA(IBAN)=DCMPLX(0.0D0,0.0D0)
 2135     CONTINUE                                                      
          DO 2230 IBAN=NBD1,NBD2                                        
            DO 1513 I=1,IIBA                                            
              I1=NBASE(I,IK)                                            
              ZTMP=ZZZ(I,IBAN,IK)*DCONJG( ZFM2( I1 ,IA ) )             
              ZTMP5=ZTMP*SSS(I,IK,KFTYPE(IA),5)
              ZTMP6=ZTMP*SSS(I,IK,KFTYPE(IA),6)
              ZTMP7=ZTMP*SSS(I,IK,KFTYPE(IA),7)
              ZTMP8=ZTMP*SSS(I,IK,KFTYPE(IA),8)
              ZTMP9=ZTMP*SSS(I,IK,KFTYPE(IA),9)
              ZTMPA=ZTMP*SSS(I,IK,KFTYPE(IA),10)
              ZFB5(IBAN)=ZFB5(IBAN)+ZTMP5
              ZFB6(IBAN)=ZFB6(IBAN)+ZTMP6
              ZFB7(IBAN)=ZFB7(IBAN)+ZTMP7
              ZFB8(IBAN)=ZFB8(IBAN)+ZTMP8
              ZFB9(IBAN)=ZFB9(IBAN)+ZTMP9
              ZFBA(IBAN)=ZFBA(IBAN)+ZTMPA
              ZFX5(IBAN)=ZFX5(IBAN)+GX(I1)*ZTMP5
              ZFX6(IBAN)=ZFX6(IBAN)+GX(I1)*ZTMP6
              ZFX7(IBAN)=ZFX7(IBAN)+GX(I1)*ZTMP7
              ZFX8(IBAN)=ZFX8(IBAN)+GX(I1)*ZTMP8
              ZFX9(IBAN)=ZFX9(IBAN)+GX(I1)*ZTMP9 
              ZFY5(IBAN)=ZFY5(IBAN)+GY(I1)*ZTMP5
              ZFY6(IBAN)=ZFY6(IBAN)+GY(I1)*ZTMP6
              ZFY7(IBAN)=ZFY7(IBAN)+GY(I1)*ZTMP7
              ZFY8(IBAN)=ZFY8(IBAN)+GY(I1)*ZTMP8
              ZFY9(IBAN)=ZFY9(IBAN)+GY(I1)*ZTMP9 
              ZFZ5(IBAN)=ZFZ5(IBAN)+GZ(I1)*ZTMP5
              ZFZ6(IBAN)=ZFZ6(IBAN)+GZ(I1)*ZTMP6
              ZFZ7(IBAN)=ZFZ7(IBAN)+GZ(I1)*ZTMP7
              ZFZ8(IBAN)=ZFZ8(IBAN)+GZ(I1)*ZTMP8
              ZFZ9(IBAN)=ZFZ9(IBAN)+GZ(I1)*ZTMP9
 1513      CONTINUE                                                    
 2230     CONTINUE                                                      
          DO 2235 IBAN=NBD1,NBD2                                        
            ZFC(IBAN,IK,IA,5)=CD*ZFB5(IBAN)
            ZFC(IBAN,IK,IA,6)=CD*ZFB6(IBAN)
            ZFC(IBAN,IK,IA,7)=CD*ZFB7(IBAN)
            ZFC(IBAN,IK,IA,8)=CD*ZFB8(IBAN)
            ZFC(IBAN,IK,IA,9)=CD*ZFB9(IBAN) 
            ZFC(IBAN,IK,IA,10)=CD*ZFBA(IBAN)
C     4/20, 1999, MODIFIED ZFX,Y,Z1 ---> ZFX,Y,Z5
            ZFORC2(IA,1)=ZFORC2(IA,1)+OCCUU(IBAN,IK)*                   
     &                   ( DCONJG(ZFX5(IBAN))*ZFC(IBAN,IK,IA,5)         
     &                    +DCONJG(ZFX6(IBAN))*ZFC(IBAN,IK,IA,6)
     &                    +DCONJG(ZFX7(IBAN))*ZFC(IBAN,IK,IA,7)
     &                    +DCONJG(ZFX8(IBAN))*ZFC(IBAN,IK,IA,8)
     &                    +DCONJG(ZFX9(IBAN))*ZFC(IBAN,IK,IA,9)  )
            ZFORC2(IA,2)=ZFORC2(IA,2)+OCCUU(IBAN,IK)*                   
     &                   ( DCONJG(ZFY5(IBAN))*ZFC(IBAN,IK,IA,5)         
     &                    +DCONJG(ZFY6(IBAN))*ZFC(IBAN,IK,IA,6)
     &                    +DCONJG(ZFY7(IBAN))*ZFC(IBAN,IK,IA,7)
     &                    +DCONJG(ZFY8(IBAN))*ZFC(IBAN,IK,IA,8)
     &                    +DCONJG(ZFY9(IBAN))*ZFC(IBAN,IK,IA,9)  ) 
            ZFORC2(IA,3)=ZFORC2(IA,3)+OCCUU(IBAN,IK)*                   
     &                   ( DCONJG(ZFZ5(IBAN))*ZFC(IBAN,IK,IA,5)         
     &                    +DCONJG(ZFZ6(IBAN))*ZFC(IBAN,IK,IA,6)
     &                    +DCONJG(ZFZ7(IBAN))*ZFC(IBAN,IK,IA,7)
     &                    +DCONJG(ZFZ8(IBAN))*ZFC(IBAN,IK,IA,8)
     &                    +DCONJG(ZFZ9(IBAN))*ZFC(IBAN,IK,IA,9)  )
 2235     CONTINUE                                                      
        END IF
C
 2110   CONTINUE                                                        
 2100 CONTINUE                                                          
!XOCL END SPREAD SUM(ZFORC2)
C!XOCL END PARALLEL
      ZCCC = 4.D0*ZI*RVOL/(2.D0*PAI)**3/FLOAT(KV3)                      
C---------------------------------------------------------              
      DO 1501 IA=1,KATM                                                 
        ZFORC2(IA,1)=ZCCC*DIMAG( ZFORC2(IA,1) )                         
        ZFORC2(IA,2)=ZCCC*DIMAG( ZFORC2(IA,2) )                         
        ZFORC2(IA,3)=ZCCC*DIMAG( ZFORC2(IA,3) )                         
 1501 CONTINUE                                                          
 3000 FORMAT(I4,6F12.8)                                                 
      RETURN                                                            
      END                                                               
C^^^^^1992 1/7 FOR STRESS ^^^^^^^^^^^                                   
      SUBROUTINE FORZFB                                                 
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                                   
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
C     VPP-PARALLEL
!XOCL SUBPROCESSOR PS(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IP=(PS,INDEX=1:KNV3,PART=BAND)
      DIMENSION SSS(KNG1,KNV3,KTYP,10)
     &         ,ZZZ(KNG1,KEG,KNV3),ZFC(KEG,KNV3,KATM,10)
!XOCL LOCAL SSS(:,/IP,:,:),ZZZ(:,:,/IP),ZFC(:,/IP,:,:)
      EQUIVALENCE (SNL,SSS),(ZAJ,ZZZ),(ZFBB,ZFC)
C                                                                       
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 1000 IA=1,KATM                                                 
      IF (NLSPD(KFTYPE(IA)).EQ.1) THEN
          LNUM = 4
      ELSE
          LNUM = 9
      END IF
      DO 1003 L=1,LNUM
      DO 1001 IK=1,KV3                                                  
      DO 1002 IBAN=NBD1,NBD2                                            
      ZFC(IBAN,IK,IA,L) =DCMPLX(0.0D0,0.0D0)                            
 1002 CONTINUE                                                          
 1001 CONTINUE                                                          
 1003 CONTINUE
 1000 CONTINUE                                                          
!XOCL END SPREAD
C!XOCL END PARALLEL
C     VPP-PARALLEL START
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 2100 IK=1,KV3                                                  
      DO 2110 IA=1,KATM                                                 
      CS=1.0D0/(WS(KFTYPE(IA))*UNIVOL)                                  
      CP=1.0D0/(WP(KFTYPE(IA))*UNIVOL)                                  
      CWL(1)=CS
      CWL(2)=CP
      CWL(3)=CP
      CWL(4)=CP
C
      IF (NLSPD(KFTYPE(IA)).EQ.1) THEN
          LNUM = 4
      ELSE
          LNUM = 9
          CD=1.0D0/(WD(KFTYPE(IA))*UNIVOL)
          CWL(5)=CD
          CWL(6)=CD
          CWL(7)=CD
          CWL(8)=CD
          CWL(9)=CD
      END IF
      DO 3222 L=1,LNUM
      DO 3200 IBAN=NBD1,NBD2                                            
C     1991 11/28  I ---> I1                                             
      DO 3510 I=1,IBA(IK)                                               
      I1  = NBASE(I,IK)                                                 
      L1  = IG1(I1)+KX1                                                 
      L2  = IG2(I1)+KY1                                                 
      L3  = IG3(I1)+KZ1                                                 
      ZTMP=ZZZ(I,IBAN,IK)*DCONJG( ZFM2( I1 ,IA ) )                     
      ZFC(IBAN,IK,IA,L)=ZFC(IBAN,IK,IA,L)+
     &                  ZTMP*SSS(I,IK,KFTYPE(IA),L)
 3510 CONTINUE                                                          
      ZFC(IBAN,IK,IA,L)=CWL(L)*ZFC(IBAN,IK,IA,L)                            
 3200 CONTINUE                                                          
 3222 CONTINUE
C
 2110 CONTINUE                                                          
 2100 CONTINUE                                                          
!XOCL END SPREAD
C!XOCL END PARALLEL
C                                                                       
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      SUBROUTINE CHAVER(IREC8,IREC9,KOPR,NOPR,KBZTYP)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
C     VPP-PARALLEL
!XOCL SUBPROCESSOR PS(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IP=(PS,INDEX=1:KNV3,PART=BAND)
      DIMENSION ZZZ(KNG1,KEG,KNV3),OCCUU(KEG,KNV3)
!XOCL LOCAL ZZZ(:,:,/IP),OCCUU(:,/IP)
      EQUIVALENCE (ZAJ,ZZZ),(OCCUP,OCCUU)
C-----ARRAYS FOR MFFT--------------------------------------------       
      DIMENSION ZC3D1(IFX2,IFY2,IFZ2)                                   
      DIMENSION IWL(8*IFX2+28),IWM(8*IFY2+28),IWN(8*IFZ2+28)            
     &         ,IWORK(2*IFX2)                                           
C================================================================       
      EQUIVALENCE (ZC1D(1),ZC3D1(1,1,1))                               
C================================================================       
C     REWIND 90                                                         
      KFT1 =  IFX2                                                      
      KFT2 =  IFY2                                                      
      KFT3 =  IFZ2                                                      
      KSUM=KFT1*KFT2*KFT3                                               
      KVOL=(KFT1-1)*KFT2*KFT3                                           
      CCC  =  2.D0*RVOL/(2.D0*PAI)**3/FLOAT(KV3)                        
      KIMG=1                                                            
C*****---- SET UP C3FFT (FAST FOURIER TRANSFORMATION) -----             
      CALL C3FFT(ZC3D1,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                
     &          ,0,0,1,IWORK,IERR)                                      
      IF (IERR.NE.0) THEN                                               
          WRITE (6,*) ' C3FFT(SET UP)]  IERR = ',IERR                   
          STOP                                                          
      END IF                                                            
      DO 7701 I=1,KSUM                                                  
        CHGB1(I)=0.0D0                                   
 7701 CONTINUE                                                          
C     VPP-PARALLEL START
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
        DO 7820 IK=1,KV3                                                
      DO 7810 IBAN=NBD1,NBD2                                            
C*******  BUILD UP CHARGE DENSITY FOR IBAN-TH BAND, IK-TH K-POINT       
C*********NAC = (IBAN-1)*KV3 + IWRT(IK)                                 
C         IF (IREC8.NE.0) NAC=IREC8*NAC-IREC8+1                         
C*********READ(80,REC=NAC) ZV1                                          
          IIBA=IBA(IK)                                                  
          DO 7830 I=1,KSUM                                              
            ZC1D(I)=DCMPLX(0.0D0,0.0D0)                                
 7830     CONTINUE                                                      
          DO 232 I=1,IIBA                                               
            I1=NBASE(I,IK)                                              
            L1=IGF1(I1)                                                 
            L2=IGF2(I1)                                                 
            L3=IGF3(I1)                                                 
            ZC3D1(L1,L2,L3) = ZZZ(I,IBAN,IK)                            
C***********ZC3D1(L1,L2,L3) = ZV1(I)                                    
  232     CONTINUE                                                      
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----                    
          CALL C3FFT(ZC3D1,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN            
     &          ,0,1,1,IWORK,IERR1)                                     
          IF (IERR1.NE.0) THEN                                          
            WRITE (6,*) ' C3FFT(IFFT C)]  IERR1 = ',IERR1               
            STOP                                                        
          END IF                                                        
C---------IREC = KV3*(IBAN-1)+IWRT(IK)                             
C            IF (IREC.LE.IRECMX) THEN                                   
C              IF (IREC9.NE.0) IREC=IREC9*IREC-IREC9+1                  
C              WRITE(90,REC=IREC)  ZC1D                                
C---------END IF                                                     
          DO 240 I=1,KSUM                                               
            CHGB1(I)=CHGB1(I)+OCCUU(IBAN,IK)*
     &               DCONJG(ZC1D(I))*ZC1D(I)                          
  240     CONTINUE                                                      
 7810 CONTINUE                                                          
 7820   CONTINUE                                                        
!XOCL END SPREAD SUM(CHGB1)
C!XOCL END PARALLEL
      DO 7900 I=1,KSUM                                                  
        ZC1D(I)=CCC*CHGB1(I)                                          
 7900 CONTINUE                                                          
C*****---- FAST FOURIER TRANSFORMATION -----                            
      CALL C3FFT(ZC3D1,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                
     &          ,0,-1,1,IWORK,IERR)                                     
      IF (IERR.NE.0) THEN                                               
          WRITE (6,*) ' C3FFT(FFT)]  IERR = ',IERR                      
          STOP                                                          
      END IF                                                            
      DENOM=1.0D0/DBLE(KVOL)                                            
      DO 30 I=1,KG                                                      
        L1=IGF1(I)                                                      
        L2=IGF2(I)                                                      
        L3=IGF3(I)                                                      
        ZCHG(I)=ZC3D1(L1,L2,L3)*DENOM                                   
   30 CONTINUE                                                          
C---------------KNGP ---> KNG KPO ---> KOPR -----------------------     
      IF (KBZTYP.GE.2) THEN                                             
          CALL CHGAVR(KOPR,NOPR,KIMG)                                                 
      END IF                                                            
C---------------------------------------------------------------        
      RETURN                                                            
      END                                                               
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7  
      SUBROUTINE CHGAVR(KOPR,NOPR,KIMG)                                                 
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      DIMENSION NPPP(KNG,KO)
      DIMENSION TAA(3,KO)      
!XOCL SUBPROCESSOR PT(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IQ=(PT,INDEX=1:KO,PART=BAND)
!XOCL LOCAL NPPP(:,/IQ),TAA(:,/IQ)
      EQUIVALENCE (NGPT,NPPP),(TAU,TAA)
      FI = 1.D0/DFLOAT(NOPR)                                            
      ZI = CMPLX(0.0D0,1.0D0)                                           
C--*-- CASE FOR REAL CHARGE                                             
          DO 200 NG = 1,KG                                              
            ZXXX(NG)    = CMPLX(0.D0,0.0D0)                            
  200     CONTINUE                                                      
!XOCL SPREAD DO /IQ
          DO 210 NO = 1,NOPR                                            
            TX = TAA(1,NO)                                              
            TY = TAA(2,NO)                                              
            TZ = TAA(3,NO)                                              
            DO 220 NG = 1,KG                                            
              NGP= NPPP(NG,NO)                                          
              FX = GX(NGP)                                              
              FY = GY(NGP)                                              
              FZ = GZ(NGP)                                              
              FP = FX*TX + FY*TY + FZ*TZ                                
              FC = EXP(ZI*FP)                                           
              ZCR= ZCHG(NGP)                                            
              ZXXX(NG)   = ZXXX(NG)    + FC*ZCR                       
  220       CONTINUE                                                    
  210     CONTINUE                                                      
!XOCL END SPREAD SUM(ZXXX)
          DO 230 NG = 1,KG                                              
            ZCHG(NG)    = ZXXX(NG)   *FI                               
  230     CONTINUE                                                      
      RETURN                                                            
      END                                                               
C     ////////////////////////////////////////////////////////////      
C     // EWALD SUMMATION FOR TOTAL ENERGY , FORCE (VEC VERSION) //      
C     //    7 OCTOBER, 1986  --->  3/15,1990                    //      
C     ////////////////////////////////////////////////////////////      
      SUBROUTINE EWVEC(ETOT2)                                           
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      PARAMETER(KEWA=KATM*(KATM+1)/2)                                   
      DIMENSION CR(KLM),CI(KLM)                                         
      DIMENSION CRX(KLM),CRY(KLM),CRZ(KLM),CLX(KLM),CLY(KLM),CLZ(KLM)   
      DIMENSION SSS(KLM),SSI(KLM)                                       
      DIMENSION ITYP(KATM)                                              
      DIMENSION FORC1(KATM,KATM,3)                                      
C
      CON=4.D0*PAI/UNIVOL                                               
      CC=2.D0/SQRT(PAI)                                                 
      ETOT2=0.D0                                                        
      IRL=KLX2*KLY2*KLZ2                                                
      IKM=KMX2*KMY2*KMZ2                                                
                MM=0                                                    
                DO 2 II=1,KTYP                                          
                I1=MM+1                                                 
                I2=MM+IATOM(II)                                         
                MM=I2                                                   
                DO 1 JJ=I1,I2                                           
C                        <------ ATTN]]]  KFTYPE <------- ]]]]]         
                ITYP(JJ)=II                                             
1               CONTINUE                                                
2               CONTINUE                                                
      IF(MM.NE.KATM)                                                    
     &  WRITE(6,*) '****WARNING "MM.NE.KATM" IN EWALD DO-LOOP 1****'    
C     ************************************************************      
C     ** THE RESULTS SHOULD BE INDEPENDENT OF THE VALUE OF ALF. **      
C     ************************************************************      
         A1X=ALTV(1,1)                                                  
         A1Y=ALTV(2,1)                                                  
         A1Z=ALTV(3,1)                                                  
                      A2X=ALTV(1,2)                                     
                      A2Y=ALTV(2,2)                                     
                      A2Z=ALTV(3,2)                                     
                                   A3X=ALTV(1,3)                        
                                   A3Y=ALTV(2,3)                        
                                   A3Z=ALTV(3,3)                        
         B1X=RLTV(1,1)                                                  
         B1Y=RLTV(2,1)                                                  
         B1Z=RLTV(3,1)                                                  
                      B2X=RLTV(1,2)                                     
                      B2Y=RLTV(2,2)                                     
                      B2Z=RLTV(3,2)                                     
                                   B3X=RLTV(1,3)                        
                                   B3Y=RLTV(2,3)                        
                                   B3Z=RLTV(3,3)                        
C                                                                       
      DO 110 J=1,KATM                                                   
      DO 111 I=1,KATM                                                   
      FORC1(I,J,1)=0.0D0                                                
      FORC1(I,J,2)=0.0D0                                                
      FORC1(I,J,3)=0.0D0                                                
  111 CONTINUE                                                          
  110 CONTINUE                                                          
C                                                                       
      MM=0                                                              
      DO 10 I1=1,KLX2                                                   
                     F1=FLOAT(I1-KLX-1)                                 
      DO 10 I2=1,KLY2                                                   
                     F2=FLOAT(I2-KLY-1)                                 
      DO 10 I3=1,KLZ2                                                   
                     F3=FLOAT(I3-KLZ-1)                                 
      MM=MM+1                                                           
      CX(MM)= F1*A1X  +  F2*A2X  +  F3*A3X                              
      CY(MM)= F1*A1Y  +  F2*A2Y  +  F3*A3Y                              
      CZ(MM)= F1*A1Z  +  F2*A2Z  +  F3*A3Z                              
   10 CONTINUE                                                          
C                                                                       
      MK=0                                                              
      DO 20 I1=1,KMX2                                                   
                     F1=FLOAT(I1-KMX-1)                                 
      DO 20 I2=1,KMY2                                                   
                     F2=FLOAT(I2-KMY-1)                                 
      DO 20 I3=1,KMZ2                                                   
                     F3=FLOAT(I3-KMZ-1)                                 
      MK=MK+1                                                           
      CKX(MK)= F1*B1X  +  F2*B2X  +  F3*B3X                            
      CKY(MK)= F1*B1Y  +  F2*B2Y  +  F3*B3Y                             
      CKZ(MK)= F1*B1Z  +  F2*B2Z  +  F3*B3Z                             
      CKR(MK)= SQRT( CKX(MK)**2 + CKY(MK)**2 + CKZ(MK)**2 )             
   20 CONTINUE                                                          
C                                                                       
      IF (IRL.NE.MM.OR.IKM.NE.MK)                                       
     &    WRITE (6,*) 'WORNING] IRL .NE. MM OR IKM .NE. MK]'            
C     STRESS 1                                                          
      DO 1000 IS=1,6                                                    
      SIGEWA(IS)=0.0D0                                                  
 1000 CONTINUE                                                          
C                                                                       
      DO 100 I=1,KATM                                                   
      DO 101 J=I,KATM                                                   
      ITE=1                                                             
      ALFF=ALF                                                         
77    CONTINUE                                                          
              AL1=SQRT(ALFF)                                            
              AL2=0.25D0/ALFF                                           
      SUM  =0.0D0                                                       
      SIG1 =0.0D0                                                       
      SIG2 =0.0D0                                                       
      SIG3 =0.0D0                                                       
      SIG4 =0.0D0                                                       
      SIG5 =0.0D0                                                       
      SIG6 =0.0D0                                                       
      SF1X =0.0D0                                                       
      SF1Y =0.0D0                                                       
      SF1Z =0.0D0                                                       
      SF1XI=0.0D0                                                       
      SF1YI=0.0D0                                                       
      SF1ZI=0.0D0                                                       
C                                                                       
      R1=CATX(I)-CATX(J)                                                
      R2=CATY(I)-CATY(J)                                                
      R3=CATZ(I)-CATZ(J)                                                
C     ///////////////////////////////////                               
C     // STEP 1;  REAL-SPACE SUMMATION //                               
C     ///////////////////////////////////                               
      DO 11 L=1,MM                                                      
      SSS(L)=0.0D0                                                      
      SSI(L)=0.0D0                                                      
      CRX(L)=CX(L)+R1                                                   
      CRY(L)=CY(L)+R2                                                   
      CRZ(L)=CZ(L)+R3                                                   
      CLX(L)=CX(L)-R1                                                   
      CLY(L)=CY(L)-R2                                                   
      CLZ(L)=CZ(L)-R3                                                   
      CR(L)=SQRT((CX(L)+R1)**2 + (CY(L)+R2)**2 + (CZ(L)+R3)**2 )        
      CI(L)=SQRT((CX(L)-R1)**2 + (CY(L)-R2)**2 + (CZ(L)-R3)**2 )        
   11 CONTINUE                                                          
C                                                                       
      DO 12 L=1,MM                                                      
      IF(CR(L).EQ.0.0D0) THEN                                           
      SUM=SUM-CC*AL1                                                    
      ELSE                                                              
      SUM=SUM+ (1.D0-ERF(AL1*CR(L)))/CR(L)                              
C-----STRESS CALCULATION                                                
          STMP = - (1.D0-ERF(AL1*CR(L)))/(CR(L)*CR(L)*CR(L))            
     &           - CC*AL1*EXP(-ALFF*CR(L)*CR(L))/(CR(L)*CR(L))          
          SIG1=SIG1+STMP*CRX(L)*CRX(L)                                  
          SIG2=SIG2+STMP*CRX(L)*CRY(L)                                  
          SIG3=SIG3+STMP*CRX(L)*CRZ(L)                                  
          SIG4=SIG4+STMP*CRY(L)*CRY(L)                                  
          SIG5=SIG5+STMP*CRY(L)*CRZ(L)                                  
          SIG6=SIG6+STMP*CRZ(L)*CRZ(L)                                  
      END IF                                                            
   12 CONTINUE                                                          
C-----FORCE CALCULATION (R & K-SPACE) --------------------------------- 
      IF (I.NE.J) THEN                                                  
      DO 14 L=1,MM                                                      
          IF (CR(L).NE.0.0D0)                                           
     &        SSS(L)=(1.0D0-ERF(AL1*CR(L)))/(CR(L)**3)+(CC*AL1          
     &                         *EXP(-ALFF*CR(L)*CR(L)))/(CR(L)*CR(L))   
          IF (CI(L).NE.0.0D0)                                           
     &        SSI(L)=(1.0D0-ERF(AL1*CI(L)))/(CI(L)**3)+(CC*AL1          
     &                         *EXP(-ALFF*CI(L)*CI(L)))/(CI(L)*CI(L))   
   14 CONTINUE                                                          
C                                                                       
      DO 15 L=1,MM                                                      
      SF1X =SF1X +CRX(L)*SSS(L)                                         
      SF1Y =SF1Y +CRY(L)*SSS(L)                                         
      SF1Z =SF1Z +CRZ(L)*SSS(L)                                         
      SF1XI=SF1XI+CLX(L)*SSI(L)                                         
      SF1YI=SF1YI+CLY(L)*SSI(L)                                         
      SF1ZI=SF1ZI+CLZ(L)*SSI(L)                                         
   15 CONTINUE                                                          
C     FORCE CALCULATION (K-SPACE) FOR OPTICAL VIBRATION                 
      DO 24 L=1,MK                                                      
      IF (CKR(L).EQ.0.0D0) GOTO 24                                      
      SSK =(CON*SIN(CKX(L)*R1+CKY(L)*R2+CKZ(L)*R3)                      
     &    *EXP(-AL2*CKR(L)*CKR(L)))/(CKR(L)*CKR(L))                     
      SF1X =SF1X +CKX(L)*SSK                                            
      SF1Y =SF1Y +CKY(L)*SSK                                            
      SF1Z =SF1Z +CKZ(L)*SSK                                            
      SF1XI=SF1XI-CKX(L)*SSK                                            
      SF1YI=SF1YI-CKY(L)*SSK                                            
      SF1ZI=SF1ZI-CKZ(L)*SSK                                            
   24 CONTINUE                                                          
C                                                                       
          FORC1(I,J,1)=                                                 
     &                ACHG(ITYP(I))*ACHG(ITYP(J))*SF1X                  
          FORC1(I,J,2)=                                                 
     &                ACHG(ITYP(I))*ACHG(ITYP(J))*SF1Y                  
          FORC1(I,J,3)=                                                 
     &                ACHG(ITYP(I))*ACHG(ITYP(J))*SF1Z                  
          FORC1(J,I,1)=                                                 
     &                ACHG(ITYP(I))*ACHG(ITYP(J))*SF1XI                 
          FORC1(J,I,2)=                                                 
     &                ACHG(ITYP(I))*ACHG(ITYP(J))*SF1YI                 
          FORC1(J,I,3)=                                                 
     &                ACHG(ITYP(I))*ACHG(ITYP(J))*SF1ZI                 
      END IF                                                            
C     /////////////////////////////////////////                         
C     // STEP 2;  RECIPROCAL-SPACE SUMMATION //                         
C     /////////////////////////////////////////                         
      DO 21 L=1,MK                                                      
C     CX = CY = CZ =0                                                   
      IF (CKR(L).EQ.0.0D0) GOTO 21                                      
      STMP=CON*EXP(-AL2*CKR(L)*CKR(L))                                  
     &   *COS(CKX(L)*R1+CKY(L)*R2+CKZ(L)*R3)/(CKR(L)*CKR(L))            
      SUM=SUM+ STMP                                                     
      SIG1=SIG1+STMP*( 2.0D0*CKX(L)*CKX(L)*(                            
     &   1.0D0/(CKR(L)*CKR(L)) + AL2 ) - 1.0D0 )                        
      SIG2=SIG2+STMP*( 2.0D0*CKX(L)*CKY(L)*(                            
     &   1.0D0/(CKR(L)*CKR(L)) + AL2 ) )                                
      SIG3=SIG3+STMP*( 2.0D0*CKX(L)*CKZ(L)*(                            
     &   1.0D0/(CKR(L)*CKR(L)) + AL2 ) )                                
      SIG4=SIG4+STMP*( 2.0D0*CKY(L)*CKY(L)*(                            
     &   1.0D0/(CKR(L)*CKR(L)) + AL2 ) - 1.0D0 )                        
      SIG5=SIG5+STMP*( 2.0D0*CKY(L)*CKZ(L)*(                            
     &   1.0D0/(CKR(L)*CKR(L)) + AL2 ) )                                
      SIG6=SIG6+STMP*( 2.0D0*CKZ(L)*CKZ(L)*(                            
     &   1.0D0/(CKR(L)*CKR(L)) + AL2 ) - 1.0D0 )                        
   21 CONTINUE                                                          
      SUM= SUM -PAI/ALFF/UNIVOL                                         
      SIG1=SIG1 +PAI/ALFF/UNIVOL                                        
      SIG4=SIG4 +PAI/ALFF/UNIVOL                                        
      SIG6=SIG6 +PAI/ALFF/UNIVOL                                        
C     ////////////                                                      
C     // STEP 3 //                                                      
C     ////////////                                                      
      IF(ITE.EQ.1) THEN                                                 
             SUM1 =SUM                                                  
             SI21 =SIG1                                                 
C------------SF2X =SF1X                                                 
C            SF2Y =SF1Y                                                 
C            SF2Z =SF1Z                                                 
C            SF2XI=SF1XI                                                
C            SF2YI=SF1YI                                                
C------------SF2ZI=SF1ZI                                                
             ALFF=2.D0*ALFF                                             
             ITE=2                                                      
C                                                                       
C            IF (MOD(ITER,ILAM).EQ.1) GOTO 77                           
             GOTO 77                                                    
      ELSE                                                              
             IF(ABS(SUM-SUM1).GT.1.D-5) WRITE(6,*)                      
     &            '******WARNING SUM.NE.SUM1***********'                
      END IF                                                            
C     WRITE (6,22) I,J,SUM1 ,SUM                                        
C     WRITE (6,22) I,J,SI21 ,SIG1                                       
C-----WRITE (6,22) I,J,SF2X ,SF1X                                       
C     WRITE (6,22) I,J,SF2Y ,SF1Y                                       
C     WRITE (6,22) I,J,SF2Z ,SF1Z                                       
C     WRITE (6,22) I,J,SF2XI,SF1XI                                      
C     WRITE (6,22) I,J,SF2YI,SF1YI                                      
C     WRITE (6,22) I,J,SF2ZI,SF1ZI                                      
      FAC=1.D0                                                          
      IF(I.EQ.J) FAC=0.5D0                                              
      ETOT2= ETOT2 +  FAC*SUM*ACHG(ITYP(I))*ACHG(ITYP(J))               
      SIGEWA(1)=SIGEWA(1)+FAC*SIG1*ACHG(ITYP(I))*ACHG(ITYP(J))          
      SIGEWA(2)=SIGEWA(2)+FAC*SIG2*ACHG(ITYP(I))*ACHG(ITYP(J))          
      SIGEWA(3)=SIGEWA(3)+FAC*SIG3*ACHG(ITYP(I))*ACHG(ITYP(J))          
      SIGEWA(4)=SIGEWA(4)+FAC*SIG4*ACHG(ITYP(I))*ACHG(ITYP(J))          
      SIGEWA(5)=SIGEWA(5)+FAC*SIG5*ACHG(ITYP(I))*ACHG(ITYP(J))          
      SIGEWA(6)=SIGEWA(6)+FAC*SIG6*ACHG(ITYP(I))*ACHG(ITYP(J))          
   22 FORMAT(1H ,2(I5,1X),'  EWALD-SUMM=',D16.8,2X,D16.8)               
  101 CONTINUE                                                          
  100 CONTINUE                                                          
      SIGEWA(1)=SIGEWA(1)/UNIVOL                                        
      SIGEWA(2)=SIGEWA(2)/UNIVOL                                        
      SIGEWA(3)=SIGEWA(3)/UNIVOL                                        
      SIGEWA(4)=SIGEWA(4)/UNIVOL                                        
      SIGEWA(5)=SIGEWA(5)/UNIVOL                                        
      SIGEWA(6)=SIGEWA(6)/UNIVOL                                        
C                                                                       
      DO 150 I=1,KATM                                                   
      FFF1(I,1)=0.0D0                                                   
      FFF1(I,2)=0.0D0                                                   
      FFF1(I,3)=0.0D0                                                   
  150 CONTINUE                                                          
C                                                                       
      DO 200 I=1,KATM                                                   
      DO 201 J=1,KATM                                                   
          FFF1(I,1)=FFF1(I,1)+FORC1(I,J,1)                              
          FFF1(I,2)=FFF1(I,2)+FORC1(I,J,2)                              
          FFF1(I,3)=FFF1(I,3)+FORC1(I,J,3)                              
  201 CONTINUE                                                          
  200 CONTINUE                                                          
C                                                                       
      IF (ITER.EQ.0) THEN                                               
      WRITE(6,*) 'TOTAL ENERGY CONTRIBUTION FROM THE  EWALD SUMM.....'  
      WRITE(6,*) ETOT2                                                  
      END IF                                                            
      RETURN                                                            
      END                                                               
C     ////////////////////////////////////////////////////////////      
C     // EWALD SUMMATION FOR TOTAL ENERGY , FORCE (VEC VERSION) //      
C     //    7 OCTOBER, 1986 ---> 3/15,1990                      //      
C     ////////////////////////////////////////////////////////////      
      SUBROUTINE EWVMD(ETOT2)                                           
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      PARAMETER(KEWA=KATM*(KATM+1)/2)                                   
      DIMENSION CR(KLM),CI(KLM)                                         
      DIMENSION CRX(KLM),CRY(KLM),CRZ(KLM),CLX(KLM),CLY(KLM),CLZ(KLM)   
      DIMENSION SSS(KLM),SSI(KLM)                                       
      DIMENSION ITYP(KATM)                                              
      DIMENSION FORC1(KATM,KATM,3)                                      
C
      CON=4.D0*PAI/UNIVOL                                               
      CC=2.D0/SQRT(PAI)                                                 
      ETOT2=0.D0                                                        
      IRL=KLX2*KLY2*KLZ2                                                
      IKM=KMX2*KMY2*KMZ2                                                
                MM=0                                                    
                DO 2 II=1,KTYP                                          
                I1=MM+1                                                 
                I2=MM+IATOM(II)                                         
                MM=I2                                                   
                DO 1 JJ=I1,I2                                           
C                        <------ ATTN]]]  KFTYPE <------- ]]]]]         
                ITYP(JJ)=II                                             
1               CONTINUE                                                
2               CONTINUE                                                
      IF(MM.NE.KATM)                                                    
     &  WRITE(6,*) '****WARNING "MM.NE.KATM" IN EWALD DO-LOOP 1****'    
C                                                                       
      DO 110 J=1,KATM                                                   
      DO 111 I=1,KATM                                                   
      FORC1(I,J,1)=0.0D0                                                
      FORC1(I,J,2)=0.0D0                                                
      FORC1(I,J,3)=0.0D0                                                
  111 CONTINUE                                                          
  110 CONTINUE                                                          
C                                                                       
      DO 100 I=1,KATM                                                   
      DO 101 J=I,KATM                                                   
      ITE=1                                                             
      ALFF=ALF                                                          
77    CONTINUE                                                          
              AL1=SQRT(ALFF)                                            
              AL2=0.25D0/ALFF                                           
      SUM  =0.0D0                                                       
      SF1X =0.0D0                                                       
      SF1Y =0.0D0                                                       
      SF1Z =0.0D0                                                       
      SF1XI=0.0D0                                                       
      SF1YI=0.0D0                                                       
      SF1ZI=0.0D0                                                       
C                                                                       
      R1=CATX(I)-CATX(J)                                                
      R2=CATY(I)-CATY(J)                                                
      R3=CATZ(I)-CATZ(J)                                                
C     ///////////////////////////////////                               
C     // STEP 1;  REAL-SPACE SUMMATION //                               
C     ///////////////////////////////////                               
      DO 11 L=1,IRL                                                     
      SSS(L)=0.0D0                                                      
      SSI(L)=0.0D0                                                      
      CRX(L)=CX(L)+R1                                                   
      CRY(L)=CY(L)+R2                                                   
      CRZ(L)=CZ(L)+R3                                                   
      CLX(L)=CX(L)-R1                                                   
      CLY(L)=CY(L)-R2                                                   
      CLZ(L)=CZ(L)-R3                                                   
      CR(L)=SQRT((CX(L)+R1)**2 + (CY(L)+R2)**2 + (CZ(L)+R3)**2 )        
      CI(L)=SQRT((CX(L)-R1)**2 + (CY(L)-R2)**2 + (CZ(L)-R3)**2 )        
   11 CONTINUE                                                          
C                                                                       
      DO 12 L=1,IRL                                                     
      IF(CR(L).EQ.0.0D0) THEN                                           
      SUM=SUM-CC*AL1                                                    
      ELSE                                                              
      SUM=SUM+ (1.D0-ERF(AL1*CR(L)))/CR(L)                              
      END IF                                                            
   12 CONTINUE                                                          
C-----FORCE CALCULATION (R-SPACE) ---------------------------------     
      IF (I.NE.J) THEN                                                  
      DO 14 L=1,IRL                                                     
          IF (CR(L).NE.0.0D0)                                           
     &        SSS(L)=(1.0D0-ERF(AL1*CR(L)))/(CR(L)**3)+(CC*AL1          
     &                         *EXP(-ALFF*CR(L)*CR(L)))/(CR(L)*CR(L))   
          IF (CI(L).NE.0.0D0)                                           
     &        SSI(L)=(1.0D0-ERF(AL1*CI(L)))/(CI(L)**3)+(CC*AL1          
     &                         *EXP(-ALFF*CI(L)*CI(L)))/(CI(L)*CI(L))   
   14 CONTINUE                                                          
C                                                                       
      DO 15 L=1,IRL                                                     
      SF1X =SF1X +CRX(L)*SSS(L)                                         
      SF1Y =SF1Y +CRY(L)*SSS(L)                                         
      SF1Z =SF1Z +CRZ(L)*SSS(L)                                         
      SF1XI=SF1XI+CLX(L)*SSI(L)                                         
      SF1YI=SF1YI+CLY(L)*SSI(L)                                         
      SF1ZI=SF1ZI+CLZ(L)*SSI(L)                                         
   15 CONTINUE                                                          
C     FORCE CALCULATION (K-SPACE) FOR OPTICAL VIBRATION                 
      DO 24 L=1,IKM                                                     
      IF (CKR(L).EQ.0.0D0) GOTO 24                                      
      SSK =(CON*SIN(CKX(L)*R1+CKY(L)*R2+CKZ(L)*R3)                      
     &    *EXP(-AL2*CKR(L)*CKR(L)))/(CKR(L)*CKR(L))                     
      SF1X =SF1X +CKX(L)*SSK                                            
      SF1Y =SF1Y +CKY(L)*SSK                                            
      SF1Z =SF1Z +CKZ(L)*SSK                                            
      SF1XI=SF1XI-CKX(L)*SSK                                            
      SF1YI=SF1YI-CKY(L)*SSK                                            
      SF1ZI=SF1ZI-CKZ(L)*SSK                                            
   24 CONTINUE                                                          
C                                                                       
          FORC1(I,J,1)=                                                 
     &                ACHG(ITYP(I))*ACHG(ITYP(J))*SF1X                  
          FORC1(I,J,2)=                                                 
     &                ACHG(ITYP(I))*ACHG(ITYP(J))*SF1Y                  
          FORC1(I,J,3)=                                                 
     &                ACHG(ITYP(I))*ACHG(ITYP(J))*SF1Z                  
          FORC1(J,I,1)=                                                 
     &                ACHG(ITYP(I))*ACHG(ITYP(J))*SF1XI                 
          FORC1(J,I,2)=                                                 
     &                ACHG(ITYP(I))*ACHG(ITYP(J))*SF1YI                 
          FORC1(J,I,3)=                                                 
     &                ACHG(ITYP(I))*ACHG(ITYP(J))*SF1ZI                 
      END IF                                                            
C     /////////////////////////////////////////                         
C     // STEP 2;  RECIPROCAL-SPACE SUMMATION //                         
C     /////////////////////////////////////////                         
      DO 21 L=1,IKM                                                     
C     CX = CY = CZ =0                                                   
      IF (CKR(L).EQ.0.0D0) GOTO 21                                      
      SUM=SUM+ CON*EXP(-AL2*CKR(L)*CKR(L))                              
     &   *COS(CKX(L)*R1+CKY(L)*R2+CKZ(L)*R3)/(CKR(L)*CKR(L))            
   21 CONTINUE                                                          
      SUM= SUM -PAI/ALFF/UNIVOL                                         
C     ////////////                                                      
C     // STEP 3 //                                                      
C     ////////////                                                      
      IF(ITE.EQ.1) THEN                                                 
             SUM1 =SUM                                                  
C------------SF2X =SF1X                                                 
C            SF2Y =SF1Y                                                 
C            SF2Z =SF1Z                                                 
C            SF2XI=SF1XI                                                
C            SF2YI=SF1YI                                                
C------------SF2ZI=SF1ZI                                                
             ALFF=2.D0*ALFF                                             
             ITE=2                                                      
      ELSE                                                              
             IF(ABS(SUM-SUM1).GT.1.D-5) WRITE(6,*)                      
     &            '******WARNING SUM.NE.SUM1***********'                
      END IF                                                            
C-----WRITE (6,22) I,J,SUM1 ,SUM                                        
C     WRITE (6,22) I,J,SF2X ,SF1X                                       
C     WRITE (6,22) I,J,SF2Y ,SF1Y                                       
C     WRITE (6,22) I,J,SF2Z ,SF1Z                                       
C     WRITE (6,22) I,J,SF2XI,SF1XI                                      
C     WRITE (6,22) I,J,SF2YI,SF1YI                                      
C     WRITE (6,22) I,J,SF2ZI,SF1ZI                                      
      FAC=1.D0                                                          
      IF(I.EQ.J) FAC=0.5D0                                              
      ETOT2= ETOT2 +  FAC*SUM*ACHG(ITYP(I))*ACHG(ITYP(J))               
   22 FORMAT(1H ,2(I5,1X),'  EWALD-SUMM=',D16.8,2X,D16.8)               
  101 CONTINUE                                                          
  100 CONTINUE                                                          
C                                                                       
      DO 150 I=1,KATM                                                   
      FFF1(I,1)=0.0D0                                                   
      FFF1(I,2)=0.0D0                                                   
      FFF1(I,3)=0.0D0                                                   
  150 CONTINUE                                                          
C                                                                       
      DO 200 I=1,KATM                                                   
      DO 201 J=1,KATM                                                   
          FFF1(I,1)=FFF1(I,1)+FORC1(I,J,1)                              
          FFF1(I,2)=FFF1(I,2)+FORC1(I,J,2)                              
          FFF1(I,3)=FFF1(I,3)+FORC1(I,J,3)                              
  201 CONTINUE                                                          
  200 CONTINUE                                                          
C                                                                       
      WRITE(6,*) 'ETOT2(EWALD) = ',ETOT2                                
      RETURN                                                            
      END                                                               
C-----------------------------------------------------------------------
      SUBROUTINE MD(IMD,KMDTP,MDFLAG,ICAR,ITERMD,ITERMX,ISTR,ISSS       
     &             ,DTUC,EPP)
C-----------------------------------------------------------------------
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
      INTEGER KMDTP(KATM)                                               
C     MOLECULAR DYNAMICS PART (SIMPLE VERLET ALGORITHM)                 
      DIMENSION VEPSI(6),EPP(3,3),VAVB(6),EPSIR(6)
      DIMENSION ALT(3,3)
C     WRITE (6,*) 'IN MD ICAR,ISTR = ',ICAR,ISTR,ISSS                   
      TKIN=0.0D0                                                        
      FCMX=0.0D0                                                        
      FDEL =-5.0D-12                                                    
      DO 10  I=1,KATM                                                   
        IF (KMDTP(I).EQ.1) THEN                                         
          FCX   = DREAL(ZFORC2(I,1))                                    
          FCY   = DREAL(ZFORC2(I,2))                                    
          FCZ   = DREAL(ZFORC2(I,3))                                    
          FCA   = SQRT(FCX**2+FCY**2+FCZ**2)                            
          IF (FCA.GT.FCMX) FCMX=FCA                                     
        END IF                                                          
  10  CONTINUE                                                          
C     WRITE (6,*) 'IN MD MDFLAG = ',MDFLAG                              
      IF(MDFLAG.EQ.0) THEN                                              
      DO 100 I=1,KATM                                                   
        IF (KMDTP(I).EQ.1) THEN                                         
          FCX   = DREAL(ZFORC2(I,1))                                    
          FCY   = DREAL(ZFORC2(I,2))                                   
          FCZ   = DREAL(ZFORC2(I,3))                                    
C**       FCA   = SQRT(FCX**2+FCY**2+FCZ**2)                            
C**       IF (FCA.GT.FCMX) FCMX=FCA                                     
          FVX=FCX*VEL(I,1)
          IF (FVX.LT.0.0D0) VEL(I,1)=0.0D0
          FVY=FCY*VEL(I,2)
          IF (FVY.LT.0.0D0) VEL(I,2)=0.0D0
          FVZ=FCZ*VEL(I,3)
          IF (FVZ.LT.0.0D0) VEL(I,3)=0.0D0
C 
          VEL(I,1) = VEL(I,1) + DTIO*FCX/AMION(I)                       
          VEL(I,2) = VEL(I,2) + DTIO*FCY/AMION(I)                       
          VEL(I,3) = VEL(I,3) + DTIO*FCZ/AMION(I)                       
          CATX(I)  = CATX(I)  + DTIO*VEL(I,1)                           
          CATY(I)  = CATY(I)  + DTIO*VEL(I,2)                           
          CATZ(I)  = CATZ(I)  + DTIO*VEL(I,3)                           
        END IF                                                          
  100 CONTINUE                                                          
C     CONSIDER BOUNDARY CONDITION                                      
C*****IF (MOD(ITER,100).EQ.0) THEN                                      
        DO 200 I=1,KATM                                                 
          P = ALINV(1,1)*CATX(I)+ALINV(1,2)*CATY(I)+ALINV(1,3)*CATZ(I)  
          Q = ALINV(2,1)*CATX(I)+ALINV(2,2)*CATY(I)+ALINV(2,3)*CATZ(I)  
          R = ALINV(3,1)*CATX(I)+ALINV(3,2)*CATY(I)+ALINV(3,3)*CATZ(I)  
C         IF (P.GE.1.0D0) P = P-1.0D0                                   
C         IF (P.LT.0.0D0) P = P+1.0D0                                   
C         IF (Q.GE.1.0D0) Q = Q-1.0D0                                   
C         IF (Q.LT.0.0D0) Q = Q+1.0D0                                   
C         IF (R.GE.1.0D0) R = R-1.0D0                                   
C         IF (R.LT.0.0D0) R = R+1.0D0                                   
C*    RECALCULATE X,Y,Z                                                 
          CATX(I)=P*ALTV(1,1)+Q*ALTV(1,2)+R*ALTV(1,3)                   
          CATY(I)=P*ALTV(2,1)+Q*ALTV(2,2)+R*ALTV(2,3)                   
          CATZ(I)=P*ALTV(3,1)+Q*ALTV(3,2)+R*ALTV(3,3)                   
          PPOS(I)=P                                                     
          QPOS(I)=Q                                                     
          RPOS(I)=R                                                     
  200   CONTINUE                                                        
C       WRITE (6,*) 'NEXT IS EPP SET]'                                  
C------------------------------------------------------------           
C     CHANGE A SHAPE OF AN UNIT CELL BY USING STRESSES      -           
C            1993 4/9                                       -           
C------------------------------------------------------------           
      IF (ICAR.EQ.2 .AND. ISTR.EQ.1) THEN                               
      WRITE (6,*) 'NOW MD SET EPP',ITER
      DO 6005 I=1,6                                                     
      VAVB(I)=0.0D0                                                     
 6005 CONTINUE                                                          
      DO 6010 IA=1,KATM                                                 
C     WRITE (6,*) 'VEL(IA) = ',IA,VEL(IA,1),VEL(IA,2),VEL(IA,3)         
      VAVB(1)=VAVB(1)+VEL(IA,1)*VEL(IA,1)/UNIVOL                        
      VAVB(2)=VAVB(2)+VEL(IA,1)*VEL(IA,2)/UNIVOL                        
      VAVB(3)=VAVB(3)+VEL(IA,1)*VEL(IA,3)/UNIVOL                        
      VAVB(4)=VAVB(4)+VEL(IA,2)*VEL(IA,2)/UNIVOL                        
      VAVB(5)=VAVB(5)+VEL(IA,2)*VEL(IA,3)/UNIVOL                        
      VAVB(6)=VAVB(6)+VEL(IA,3)*VEL(IA,3)/UNIVOL                        
 6010 CONTINUE                                                          
      IF (ISSS.EQ.0) THEN                                               
          DO 6000 I=1,6                                                 
C         VEPSI(I)=-DTUC*(SIGSTR(I)+VAVB(I)-PEX(I))                     
          EPSIR(I)=-DTUC*(SIGSTR(I)+VAVB(I)-PEX(I))                     
C         EPSIR(I)= DTUC*VEPSI(I)                                       
 6000 CONTINUE                                                          
          ISSS=1                                                        
      END IF                                                            
      DO 6100 I=1,6                                                     
C     CHKVF=VEPSI(I)*(SIGSTR(I)+VAVB(I)-PEX(I))                         
      CHKVF=EPSIR(I)*(SIGSTR(I)+VAVB(I)-PEX(I))                         
        IF (CHKVF.GE.0.0D0)  THEN                                       
            EPSIR(I)=0.0D0                                              
C            VEPSI(I)=0.0D0                                              
        END IF                                                          
C           VEPSI(I)=VEPSI(I) -DTUC*(SIGSTR(I)+VAVB(I)-PEX(I))          
            EPSIR(I)=         -DTUC*(SIGSTR(I)+VAVB(I)-PEX(I))          
C           EPSIR(I)=EPSIR(I) +DTUC*VEPSI(I)                            
 6100 CONTINUE                                                          
        EPP(1,1)=EPSIR(1) + 1.0D0                                       
        EPP(1,2)=EPSIR(2)
        EPP(1,3)=EPSIR(3)
        EPP(2,1)=EPSIR(2)
        EPP(2,2)=EPSIR(4) + 1.0D0                                       
        EPP(2,3)=EPSIR(5)
        EPP(3,1)=EPSIR(3)
        EPP(3,2)=EPSIR(5)
        EPP(3,3)=EPSIR(6) + 1.0D0                                       
        WRITE (6,*) 'EPP(1,1) = ',EPP(1,1),VAVB(1)                      
        WRITE (6,*) 'EPP(1,2) = ',EPP(1,2),VAVB(2)                      
        WRITE (6,*) 'EPP(1,3) = ',EPP(1,3),VAVB(3)                      
        WRITE (6,*) 'EPP(2,1) = ',EPP(2,1),VAVB(2)                      
        WRITE (6,*) 'EPP(2,2) = ',EPP(2,2),VAVB(4)                      
        WRITE (6,*) 'EPP(2,3) = ',EPP(2,3),VAVB(5)                      
        WRITE (6,*) 'EPP(3,1) = ',EPP(3,1),VAVB(3)                      
        WRITE (6,*) 'EPP(3,2) = ',EPP(3,2),VAVB(5)                      
        WRITE (6,*) 'EPP(3,3) = ',EPP(3,3),VAVB(6)                      
 5000 CONTINUE                                                          
C     CHANGE AXES OF A UNIT CELL                                        
      DO 6200 I=1,3                                                     
      ALT(I,1)=EPP(I,1)*ALTV(1,1)+EPP(I,2)*ALTV(2,1) 
     &        +EPP(I,3)*ALTV(3,1)                   
      ALT(I,2)=EPP(I,1)*ALTV(1,2)+EPP(I,2)*ALTV(2,2)
     &        +EPP(I,3)*ALTV(3,2)                   
      ALT(I,3)=EPP(I,1)*ALTV(1,3)+EPP(I,2)*ALTV(2,3)
     &        +EPP(I,3)*ALTV(3,3)                   
 6200 CONTINUE                                                          
C
      DO 6210 I=1,3
      ALTV(I,1)=ALT(I,1)
      ALTV(I,2)=ALT(I,2)
      ALTV(I,3)=ALT(I,3)
 6210 CONTINUE
C                                                                       
      DO 6500 IA=1,KATM
      P=PPOS(IA)
      Q=QPOS(IA)
      R=RPOS(IA)
      CATX(IA)=P*ALTV(1,1)+Q*ALTV(1,2)+R*ALTV(1,3)
      CATY(IA)=P*ALTV(2,1)+Q*ALTV(2,2)+R*ALTV(2,3)
      CATZ(IA)=P*ALTV(3,1)+Q*ALTV(3,2)+R*ALTV(3,3) 
 6500 CONTINUE 
C
      END IF                                                            
C                                                                       
C      IF(MOD(ITER,10).EQ.0) THEN                                       
        WRITE (18,*) 'ITER = ',ITER                                     
        WRITE (18,300)(I,CATX(I),CATY(I),CATZ(I),                       
     &                   PPOS(I),QPOS(I),RPOS(I)  ,I=1,KATM)            
C      END IF                                                           
      IF((FCMX.LT.FDEL*0.5D0).OR.(ITER.GE.ITERMX)) THEN                 
        MDFLAG  = 1                                                     
        ICAR    = 0                                                     
        ITERMD  = ITER                                                  
                WRITE(6,*) ' ITERMD= ',ITERMD                           
      END IF                                                            
      END IF                                                            
      IF((FCMX.GT.FDEL).AND.(ITER.LE.ITERMX)) THEN                      
        MDFLAG  = 0                                                     
        ICAR    = 2                                                     
        ITERMD  = 5000                                                  
      END IF                                                            
  300 FORMAT(' ',I3,3F15.10,3F10.6)                                     
  302 FORMAT(' ',I3,3F20.10)                                            
      WRITE(21,305) ITER,FCMX,MDFLAG                                    
  305 FORMAT(' ',I5,F12.6,3X,I3,'ITER AND FCMX')                        
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                        
      SUBROUTINE MSD(IREC8,IREC9,IMD,ICAR,ISTR)                         
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                        
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
C     VPP-PARALLEL
!XOCL SUBPROCESSOR PS(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IP=(PS,INDEX=1:KNV3,PART=BAND)
      DIMENSION SSS(KNG1,KNV3,KTYP,10),ZZZ(KNG1,KEG,KNV3)
     &         ,ZFC(KEG,KNV3,KATM,10)
      DIMENSION EKK(KEG,KNV3),OCCUU(KEG,KNV3),RAA(KNG1,KNV3)
!XOCL LOCAL SSS(:,/IP,:,:),ZZZ(:,:,/IP),ZFC(:,/IP,:,:)
!XOCL LOCAL EKK(:,/IP),OCCUU(:,/IP),RAA(:,/IP)
      EQUIVALENCE (SNL,SSS),(ZAJ,ZZZ),(ZFBB,ZFC)
     &           ,(EKO,EKK),(OCCUP,OCCUU),(RAK,RAA)
C     EIGEN-VALUE PROBLEM                                               
      DIMENSION ZDEVC(KNG1,KEG)                                         
     &         ,EMAS(KEG)
C-----ARRAYS FOR MFFT--------------------------------------------       
      DIMENSION ZEV3C(IFX2,IFY2,IFZ2),ZV3D(IFX2,IFY2,IFZ2)              
      DIMENSION IWL(8*IFX2+28),IWM(8*IFY2+28),IWN(8*IFZ2+28)            
     &         ,IWORK(2*IFX2)                                           
C================================================================       
      EQUIVALENCE (ZV1D(1),ZV3D(1,1,1)),(ZC1D(1),ZEV3C(1,1,1))          
C==== ATTN KX1 OR KX11 ==========================================       
C     REWIND 90                                                         
      KFT1 =  IFX2                                                      
      KFT2 =  IFY2                                                      
      KFT3 =  IFZ2                                                      
      KSUM=KFT1*KFT2*KFT3                                               
      KVOL=(KFT1-1)*KFT2*KFT3                                           
      IF (ITER.EQ.2) WRITE (6,*) 'DTIME = ',DTIM                        
C---------- MASS OF ELECTRON -----                                      
      DO 4 I=1,KG                                                       
        ZCHGO(I) = ZCHG(I)                                              
        ZCHG(I)  = DCMPLX(0.0D0,0.0D0)                                  
    4 CONTINUE                                                          
C     FOR STRESS (CALL KBINT & FORZFB 1992 1/7)                         
      IF(ISTR.EQ.1) THEN                                
      WRITE (6,*) 'KBINT ITER>IMD IN MSD'                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                               
      CALL KBINT
C^^^^^    STRESS     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                
      CALL FORZFB                                                       
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                               
      END IF                                                            
C     ////////////////////////////////////                              
C     // CALCULATE THE VERLET ALGORITHM //                              
C     ////////////////////////////////////                              
      ZVP(1) = ZVXC(1) + ZPSCC(1)                                       
      DO 2001 I=2,KG                                                    
        ZVP(I) = ZVXC(I) + ZPSCC(I) + PAI4*ZCHGO(I)*RGG(I)
 2001 CONTINUE                                                          
C*****---- SET UP C3FFT (FAST FOURIER TRANSFORMATION) -----             
      CALL C3FFT(ZEV3C,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                
     &          ,0,0,1,IWORK,IERR)                                      
      IF (IERR.NE.0) THEN                                               
          WRITE (6,*) ' C3FFT(SET UP)]  IERR = ',IERR                   
          STOP                                                          
      END IF                                                            
      DO 234 I=1,KSUM                                                   
        ZV1D(I)=DCMPLX(0.0D0,0.0D0)                                     
  234 CONTINUE                                                          
      DO 233 I=1,KG                                                     
        LF1=IGF1(I)                                                     
        LF2=IGF2(I)                                                     
        LF3=IGF3(I)                                                     
        ZV3D(LF1,LF2,LF3) = ZVP(I)                                      
  233 CONTINUE                                                          
      CALL TIME(ITIME1,'      ',1)                                      
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----                    
      CALL C3FFT(ZV3D ,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                
     &          ,0,1,1,IWORK,IERR2)                                     
      CALL TIME(ITIME1,'C3FFT ',2)                                      
      IF (IERR2.NE.0) THEN                                              
          WRITE (6,*) ' C3FFT(IFFT V)]  IERR2 = ',IERR2                 
          STOP                                                          
      END IF                                                            
C     WRITE(6,*) 'ZV1D = '                                              
C     DO 27 I=1,20                                                      
C       WRITE(6,*) ZV1D(I)                                              
C  27 CONTINUE                                                          
C-----K-POINT LOOP---------------------------------------------         
C     VPP-PARALLEL START 2
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 100 NNN=1,KV3                                                  
C                                  IWRT(NNN) =NNN                       
                                   AKX = VX(NNN)                        
                                   AKY = VY(NNN)                        
                                   AKZ = VZ(NNN)                        
                                  IIBA = IBA(NNN)                       
      DO 3 IBAN=NBD1,NBD2                                               
      EMAS(IBAN)=0.0D0                                                  
    3 CONTINUE                                                          
      DO 2002 J=1,IIBA                                                  
        X(J) = 0.0D0
        VKIN =                                                       
     &           ( (AKX+GX(NBASE(J,NNN)))**2                            
     &           + (AKY+GY(NBASE(J,NNN)))**2                            
     &           + (AKZ+GZ(NBASE(J,NNN)))**2 )/2.D0                     
        Y3(J)   =  VKIN + DREAL(ZVP(1))
 2002 CONTINUE                                                          
      DO 7000 ITY=1,KTYP                                                
        CS=1.0D0/(WS(ITY)*UNIVOL)                                       
        CP=1.0D0/(WP(ITY)*UNIVOL)                                       
        CWL(1)=CS
        CWL(2)=CP
        CWL(3)=CP
        CWL(4)=CP
        IF (NLSPD(ITY).EQ.1) THEN
            LNUM = 4
        ELSE
            LNUM = 9
            CD=1.0D0/(WD(ITY)*UNIVOL)
            CWL(5)=CD
            CWL(6)=CD
            CWL(7)=CD
            CWL(8)=CD
            CWL(9)=CD
        END IF
C
        DO 7111 L=1,LNUM
        DO 7010 I=1,IIBA                                                
          TMPP       = CWL(L)*SSS(I,NNN,ITY,L)**2
          Y3(I)  = Y3(I) + TMPP*DFLOAT(IATOM(ITY))         
          X(I)   = X(I)  + TMPP*DFLOAT(IATOM(ITY))         
 7010   CONTINUE                                                        
 7111   CONTINUE
 7000 CONTINUE                                                          
C-----BAND LOOP----------------------------------------------           
      DO 200 IBAN=NBD1,NBD2                                             
C*********IREC = KV3*(IBAN-1)+IWRT(NNN)                                 
C         IF (IREC8.NE.0) IREC=IREC8*IREC-IREC8+1                       
C         READ(80,REC=IREC) ZV1                                         
C         IREC = KV3*(IBAN-1)+IWRT(NNN)                                 
C         IF (IREC.LE.1000) THEN
C            IF (IREC.LE.IRECMX) THEN                                   
C         IF (IREC9.NE.0) IREC=IREC9*IREC-IREC9+1                       
C           READ(90,REC=IREC) ZRC
C        ELSE                                                           
C        READ(90,REC=IREC) ZRC                                         
C---------------------------------------90.1.22 Y.M.----                
      DO 3231 I=1,KSUM                                                  
        ZC1D(I)=CMPLX(0.0D0,0.0D0)                                      
 3231 CONTINUE                                                          
      DO 3232 I=1,IIBA                                                  
        I1=NBASE(I,NNN)                                                 
        L1=IGF1(I1)                                                     
        L2=IGF2(I1)                                                     
        L3=IGF3(I1)                                                     
        ZEV3C(L1,L2,L3) = ZZZ(I,IBAN,NNN)                               
C*******ZEV3C(L1,L2,L3) = ZV1(I)                                        
 3232 CONTINUE                                                          
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----                    
      CALL C3FFT(ZEV3C,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                
     &          ,0,1,1,IWORK,IERR1)                                     
      IF (IERR1.NE.0) THEN                                              
          WRITE (6,*) ' C3FFT(IFFT C)]  IERR1 = ',IERR1                 
          STOP                                                          
      END IF                                                            
CXX   DO 3240 I=1,KSUM                                                  
CXX     ZRC(I)=ZC1D(I)                                                  
C3240 CONTINUE                                                          
C        END IF                                                         
C     ----- CALCULATE THE MATRIX ELEMENTS -----                         
      DO 71 I=1,KNG1
C        ZPNL(I) = - ZZZ(I,IBAN,NNN)*X(I)
        Y1(I) = - DREAL(ZZZ(I,IBAN,NNN)*X(I))
        Y2(I) = - DIMAG(ZZZ(I,IBAN,NNN)*X(I))
   71 CONTINUE                                                          
      DO 401 IA=1,KATM                                                  
C
      IF (NLSPD(KFTYPE(IA)).EQ.1) THEN
          LNUM = 4
      ELSE
          LNUM = 9
      END IF
C
        DO 411 L =1,LNUM
        DO 400 I =1,IIBA                                                
          ZTMP = SSS(I,NNN,KFTYPE(IA),L)                                
     &         * ZFC(IBAN,NNN,IA,L)                                     
C          ZPNL(I)   = ZPNL(I)  + ZTMP*ZFM2(NBASE(I,NNN),IA)             
          Y1(I) = Y1(I) + DREAL(ZTMP*ZFM2(NBASE(I,NNN),IA))
          Y2(I) = Y2(I) + DIMAG(ZTMP*ZFM2(NBASE(I,NNN),IA))
  400   CONTINUE                                                        
  411   CONTINUE
  401 CONTINUE                                                          
      DO 240 I=1,KSUM                                                   
C       ZC1D(I)=ZV1D(I)*ZRC(I)                                          
        ZC1D(I)=ZV1D(I)*ZC1D(I)                                         
  240 CONTINUE                                                          
C*****---- FAST FOURIER TRANSFORMATION -----                            
      CALL C3FFT(ZEV3C,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN                
     &          ,0,-1,1,IWORK,IERR)                                     
      IF (IERR.NE.0) THEN                                               
          WRITE (6,*) ' C3FFT(FFT)]  IERR = ',IERR                      
          STOP                                                          
      END IF                                                            
      DENOM=1.0D0/DBLE(KVOL)                                            
      DO 30 I=1,IIBA                                                    
          I1=NBASE(I,NNN)                                               
          L1=IGF1(I1)                                                   
          L2=IGF2(I1)                                                   
          L3=IGF3(I1)                                                   
          ZPGG=ZEV3C(L1,L2,L3)*DENOM                                 
      WDI  =(Y3(I)-EKK(IBAN,NNN))
C      ZROF =(ZPGG+ZPNL(I)-ZVP(1)*ZZZ(I,IBAN,NNN))
      ZROF =(ZPGG+Y1(I)+ZI*Y2(I)-ZVP(1)*ZZZ(I,IBAN,NNN))
      ZDEVC(I,IBAN)=-(WDI*ZZZ(I,IBAN,NNN)+ZROF)
      EMAS(IBAN)=EMAS(IBAN)+DCONJG(ZZZ(I,IBAN,NNN))*ZDEVC(I,IBAN)
      ESHI=DEXP(-WDI*DTIM)
      ZZZ(I,IBAN,NNN)=ZZZ(I,IBAN,NNN)*ESHI+(ESHI-1.0D0)*ZROF/WDI  
   30 CONTINUE                                                          
C-------------------------------------------------------------------    
  200 CONTINUE                                                          
C     >>>>> AFTER C(T+DT)=C(T)+DC(T) <<<<<                              
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                               
C     GRAM-SCHMIDT ORTHOGONALIZATION                                    
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                               
C      IF (NNN.EQ.1) CALL TIME(ITIME1,'      ',1)                       
C      CALL GRAMMD(NNN,IIBA)                                            
C      IF (NNN.EQ.1) CALL TIME(ITIME1,'GRAM  ',2)                       
      DO 1000 IBAN=NBD1,NBD2-1                                          
        ZNRM=DCMPLX(0.0D0,0.0D0)                                        
        DO 1010 I=1,IIBA                                                
          ZNRM=ZNRM+DCONJG(ZZZ(I,IBAN,NNN))*ZZZ(I,IBAN,NNN)             
 1010   CONTINUE                                                        
        DENOM=1.0D0/DSQRT(DREAL(ZNRM))                                  
        DO 1020 I=1,IIBA                                                
          ZZZ(I,IBAN,NNN)=ZZZ(I,IBAN,NNN)*DENOM                         
 1020   CONTINUE                                                        
        DO 1030 JBAN=IBAN+1,NBD2                                        
          ZNRM=DCMPLX(0.0D0,0.0D0)                                      
          DO 1040 I=1,IIBA                                              
            ZNRM=ZNRM+DCONJG(ZZZ(I,IBAN,NNN))*ZZZ(I,JBAN,NNN)           
 1040     CONTINUE                                                      
          DO 1050 I=1,IIBA                                              
            ZZZ(I,JBAN,NNN)=ZZZ(I,JBAN,NNN)-ZNRM*ZZZ(I,IBAN,NNN)        
 1050     CONTINUE                                                      
 1030   CONTINUE                                                        
 1000 CONTINUE                                                          
      ZNRM=DCMPLX(0.0D0,0.0D0)                                         
      DO 1100 I=1,IIBA                                                  
        ZNRM=ZNRM+DCONJG(ZZZ(I,NBD2,NNN))*ZZZ(I,NBD2,NNN)               
 1100 CONTINUE                                                          
      DENOM=1.0D0/DSQRT(DREAL(ZNRM))                                    
      DO 1110 I=1,IIBA                                                  
        ZZZ(I,NBD2,NNN)=ZZZ(I,NBD2,NNN)*DENOM                           
 1110 CONTINUE                                                          
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                               
C     CHECK                                                             
C**************************************************************         
C     IF (ITER.GT.40) THEN                                              
C     DO 700 IBAN=NBD1,NBD2                                             
C       DO 720 JBAN=NBD1,IBAN                                           
C         ZNORM=DCMPLX(0.0D0,0.0D0)                                     
C         DO 7001 I=1,IIBA                                              
C           ZNORM=ZNORM+DCONJG(ZAJ(I,IBAN,NNN))*ZAJ(I,JBAN,NNN)         
C7001     CONTINUE                                                      
C         RNOR=DSQRT(CDABS(ZNORM))                                      
C         IF (IBAN.EQ.JBAN) THEN                                        
C             IF(DABS(1.0D0-RNOR).GT.1.0D-7) THEN                       
C               WRITE(6,*) 'NRM.NE.1.0 I= ',IBAN,'NRM=',RNOR            
C             END IF                                                    
C           ELSE                                                        
C             IF(RNOR.GT.1.0D-7) THEN                                   
C               WRITE(6,*) 'NRM.NE.0.0 I J= ',IBAN,JBAN                 
C    &                           ,'NRM= ',RNOR                          
C             END IF                                                    
C             IFLAG=1                                                   
C         END IF                                                        
C 720   CONTINUE                                                        
C 700 CONTINUE                                                          
C     END IF                                                            
C******************************************************************     
      DO 261 IBAN=NBD1,NBD2                                             
      DO 260 I=1,IIBA                                                   
      EKK(IBAN,NNN)=EKK(IBAN,NNN)                                       
     &             -2.0D0*DCONJG(ZZZ(I,IBAN,NNN))*ZDEVC(I,IBAN)         
  260 CONTINUE                                                          
      EKK(IBAN,NNN)=EKK(IBAN,NNN)+EMAS(IBAN)
C      EG(IBAN)=EKK(IBAN,NNN)                                           
  261 CONTINUE                                                          
      DO 262 IBAN=NBD1,NBD2-1                                           
      DO 262 JBAN=IBAN+1,NBD2                                           
         IF(EKK(JBAN,NNN).LT.EKK(IBAN,NNN)) THEN                        
            EE =EKK(IBAN,NNN)                                           
            EKK(IBAN,NNN)=EKK(JBAN,NNN)                                 
            EKK(JBAN,NNN)=EE                                            
            DO 270 IV=1,IIBA                                            
            ZTV              = ZZZ(IV,IBAN,NNN)
            ZZZ(IV,IBAN,NNN) = ZZZ(IV,JBAN,NNN)                         
            ZZZ(IV,JBAN,NNN) = ZTV                                      
  270 CONTINUE                                                          
         END IF                                                         
  262 CONTINUE                                                          
C      DO 263 IBAN=NBD1,NBD2                                            
C      EKK(IBAN,NNN)=EG(IBAN)                                           
C     WRITE(6,*) 'NNN = ',NNN,'IBAN = ',IBAN                            
C     WRITE(6,*) EKO(IBAN,NNN)                                          
C  263 CONTINUE                                                         
C     OUTPUT EIGENVECTORS ON FILE 80 (DIRECT-ACCESS) (KG1 <---> IIBA)   
      DO 231 IBAN = NBD1,NBD2                                           
C**************IREC = KV3*(IBAN-1)+IWRT(NNN)                            
C              IF (IREC8.NE.0) IREC=IREC8*IREC-IREC8+1                  
C**************WRITE(80,REC=IREC)  ZV1                                  
  231 CONTINUE                                                          
  100 CONTINUE                                                          
!XOCL END SPREAD
C!XOCL END PARALLEL
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                 
      SUBROUTINE STRESS(IPCC,SCHGPC,ETOT1,TOTCH,EPC,PCM    
     &                     ,KOPR,NBZTYP)
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                 
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
C     VPP-PARALLEL
!XOCL SUBPROCESSOR PS(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IP=(PS,INDEX=1:KNV3,PART=BAND)
      DIMENSION ZZZ(KNG1,KEG,KNV3),OCCUU(KEG,KNV3)
!XOCL LOCAL ZZZ(:,:,/IP),OCCUU(:,/IP)
      EQUIVALENCE (ZAJ,ZZZ),(OCCUP,OCCUU)
      DIMENSION OO(3,3,KO)
!XOCL INDEX PARTITION IQ=(PS,INDEX=1:KO,PART=BAND)
!XOCL LOCAL OO(:,:,/IQ)
      EQUIVALENCE (OP,OO)
C
      DIMENSION SSS(3,3),TTT(3,3)                                       
C     STRESS CALCULATION FOR KINETIC, HARTREE, EX-COR, LOCAL AND ETOT1  
      CCC =2.D0*RVOL/(2.D0*PAI)**3/FLOAT(KV3)                           
      DO 1000 IS=1,6                                                    
      SIGSTR(IS)=0.0D0                                                  
 1000 CONTINUE                                                          
      RGG(1)=0.0D0                                                       
      ZXXX(1)=DCMPLX(0.0D0,0.0D0)
      ZYYY(1)=DCMPLX(0.0D0,0.0D0)
      ZZZZ(1)=DCMPLX(0.0D0,0.0D0)
      DO 1001 I=2,KG                                                    
      ZXXX(I)=DCMPLX(0.0D0,0.0D0)
      ZYYY(I)=DCMPLX(0.0D0,0.0D0)
      ZZZZ(I)=DCMPLX(0.0D0,0.0D0)
      RGG(I)=1.0D0/(GR(I)*GR(I))                                         
 1001 CONTINUE                                                          
C     KINETIC PART                                                      
C     VPP-PARALLEL START
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 1100 IK=1,KV3                                                  
C                                  IWRT(IK) =IK                         
                                   AKX = VX(IK)                         
                                   AKY = VY(IK)                         
                                   AKZ = VZ(IK)                         
      DO 1110 IBAN=NBD1,NBD2                                            
      CW=CCC*OCCUU(IBAN,IK)                                             
      DO 1120 I=1,IBA(IK)                                               
      I1  = NBASE(I,IK)                                                 
      STMP=CW*DCONJG(ZZZ(I,IBAN,IK))*ZZZ(I,IBAN,IK)                     
      SIGSTR(1)=SIGSTR(1)-STMP*(AKX+GX(I1))*(AKX+GX(I1))                
      SIGSTR(2)=SIGSTR(2)-STMP*(AKX+GX(I1))*(AKY+GY(I1))                
      SIGSTR(3)=SIGSTR(3)-STMP*(AKX+GX(I1))*(AKZ+GZ(I1))                
      SIGSTR(4)=SIGSTR(4)-STMP*(AKY+GY(I1))*(AKY+GY(I1))                
      SIGSTR(5)=SIGSTR(5)-STMP*(AKY+GY(I1))*(AKZ+GZ(I1))                
      SIGSTR(6)=SIGSTR(6)-STMP*(AKZ+GZ(I1))*(AKZ+GZ(I1))                
 1120 CONTINUE                                                          
 1110 CONTINUE                                                          
 1100 CONTINUE                                                          
!XOCL END SPREAD SUM(SIGSTR)
C!XOCL END PARALLEL
      SSS1=SIGSTR(1)                                                    
      SSS2=SIGSTR(2)                                                    
      SSS3=SIGSTR(3)                                                    
      SSS4=SIGSTR(4)                                                    
      SSS5=SIGSTR(5)                                                    
      SSS6=SIGSTR(6)                                                    
      WRITE (6,*) 'TOTAL STR 1  1 = ',SIGSTR(1)                         
      WRITE (6,*) '             2 = ',SIGSTR(2)                         
      WRITE (6,*) '             3 = ',SIGSTR(3)                         
      WRITE (6,*) '             4 = ',SIGSTR(4)                         
      WRITE (6,*) '             5 = ',SIGSTR(5)                         
      WRITE (6,*) '             6 = ',SIGSTR(6)                         
C     HARTREE PART   I=1 -----> GR(1)=0.0D0 : RGG(1)=INFINITE            
C     EXCHANGE-CORRELATION + LOCAL PS PARTS                             
      DO 1250 I=1,KG                                                    
      ZVP(I)=ZVXC(I)                                                  
 1250 CONTINUE                                                          
      CALL XCFFT(1,SCHGPC)                                 
      CALL XCFFT(2,SCHGPC)                                 
      IF (IPCC.EQ.1) CALL XSTPC(1,SCHGPC)                   
      STMQ=DCONJG(ZEXC(1))*(ZCHG(1)+ZRHPC(1))-DCONJG(ZVXC(1))*ZCHG(1)   
     &    -DCONJG(ZVXC(1))*ZRHPC(1)                                     
      SS1=-DCONJG(ZEXC(1))*ZRRPC(1)                                     
      SSX=-DCONJG(ZXXX(1))*(ZCHG(1)+ZRHPC(1))                           
      SSY=-DCONJG(ZYYY(1))*(ZCHG(1)+ZRHPC(1))                           
      SSZ=-DCONJG(ZZZZ(1))*(ZCHG(1)+ZRHPC(1))                           
      SIGSTR(1)=SIGSTR(1)+STMQ+SS1*GX(1)*GX(1)+SSX                      
      SIGSTR(4)=SIGSTR(4)+STMQ+SS1*GY(1)*GY(1)+SSY                      
      SIGSTR(6)=SIGSTR(6)+STMQ+SS1*GZ(1)*GZ(I)+SSZ                      
      TTT1=SIGSTR(1)-SSS1                                               
      TTT2=SIGSTR(2)-SSS2                                               
      TTT3=SIGSTR(3)-SSS3                                               
      TTT4=SIGSTR(4)-SSS4                                               
      TTT5=SIGSTR(5)-SSS5                                               
      TTT6=SIGSTR(6)-SSS6                                               
      SSS1=SIGSTR(1)                                                    
      SSS2=SIGSTR(2)                                                    
      SSS3=SIGSTR(3)                                                    
      SSS4=SIGSTR(4)                                                    
      SSS5=SIGSTR(5)                                                    
      SSS6=SIGSTR(6)                                                    
      WRITE (6,*) 'TOTAL STR 2  1 = ',TTT1                              
      WRITE (6,*) '             2 = ',TTT2                              
      WRITE (6,*) '             3 = ',TTT3                              
      WRITE (6,*) '             4 = ',TTT4                              
      WRITE (6,*) '             5 = ',TTT5                              
      WRITE (6,*) '             6 = ',TTT6                              
      STX=0.0D0                                                         
      STY=0.0D0                                                         
      STZ=0.0D0                                                         
      STT=0.0D0                                                         
      STU=0.0D0
      STV=0.0D0
      STW=0.0D0
      DO 1200 I=2,KG                                                    
      STMP=0.5D0*PAI4*DCONJG(ZCHG(I))*ZCHG(I)*RGG(I)                     
      STMQ=DCONJG(ZEXC(I))*(ZCHG(I)+ZRHPC(I))-DCONJG(ZVXC(I))*ZCHG(I)   
     &    -DCONJG(ZVXC(I))*ZRHPC(I)                                     
      STM1= DCONJG(ZEXC(I))*ZRRPC(I)                                    
      STMX= DCONJG(ZXXX(I))*(ZCHG(I)+ZRHPC(I))                          
      STMY= DCONJG(ZYYY(I))*(ZCHG(I)+ZRHPC(I))                          
      STMZ= DCONJG(ZZZZ(I))*(ZCHG(I)+ZRHPC(I))                          
C      STX=STX+DREAL(DCONJG(ZVXC(I))*ZRRPC(I))*GX(I)*GX(I)              
C      STY=STY+DREAL(DCONJG(ZEXC(I))*ZRRPC(I))*GX(I)*GX(I)              
C      STZ=STZ+STMX                                                     
C      STT=STT+STMY                                                     
      STX=STX+STMP*(2.0D0*GX(I)*GX(I)*RGG(I)-1.0D0)+STMQ     
     &   -STMS*2.0D0*GX(I)*GX(I)-STMR-STM1*GX(I)*GX(I)-STMX       
      STY=STY+STMP*(2.0D0*GX(I)*GX(I)*RGG(I)-1.0D0)
      STZ=STZ+STMQ
      STT=STT-STMS*2.0D0*GX(I)*GX(I)
      STU=STU-STMR
      STV=STV-STM1*GX(I)*GX(I)
      STW=STW-STMX
      STMR=ZPSCC(I)*DCONJG(ZCHG(I))                                     
      STMS=ZDSCC(I)*DCONJG(ZCHG(I))                                     
      SIGSTR(1)=SIGSTR(1)+STMP*(2.0D0*GX(I)*GX(I)*RGG(I)-1.0D0)+STMQ     
     &         -STMS*2.0D0*GX(I)*GX(I)-STMR-STM1*GX(I)*GX(I)-STMX       
      SIGSTR(2)=SIGSTR(2)+STMP*(2.0D0*GX(I)*GY(I)*RGG(I))                
     &         -STMS*2.0D0*GX(I)*GY(I)                                  
      SIGSTR(3)=SIGSTR(3)+STMP*(2.0D0*GX(I)*GZ(I)*RGG(I))                
     &         -STMS*2.0D0*GX(I)*GZ(I)                                  
      SIGSTR(4)=SIGSTR(4)+STMP*(2.0D0*GY(I)*GY(I)*RGG(I)-1.0D0)+STMQ     
     &         -STMS*2.0D0*GY(I)*GY(I)-STMR-STM1*GY(I)*GY(I)-STMY       
      SIGSTR(5)=SIGSTR(5)+STMP*(2.0D0*GY(I)*GZ(I)*RGG(I))                
     &         -STMS*2.0D0*GY(I)*GZ(I)                                  
      SIGSTR(6)=SIGSTR(6)+STMP*(2.0D0*GZ(I)*GZ(I)*RGG(I)-1.0D0)+STMQ     
     &         -STMS*2.0D0*GZ(I)*GZ(I)-STMR-STM1*GZ(I)*GZ(I)-STMZ       
 1200 CONTINUE                                                          
      WRITE (6,*) 'ZVX,ZEX = ',STX,STY                                  
      WRITE (6,*) 'STM1(XY)= ',STZ,STT                                  
      WRITE (6,*) 'STU,V,W = ',STU,STV,STW
      STOT=STY+STZ+STT+STU+STV+STW
      WRITE (6,*) 'STOT    = ',STOT
      TTT1=SIGSTR(1)-SSS1                                               
      TTT2=SIGSTR(2)-SSS2                                               
      TTT3=SIGSTR(3)-SSS3                                               
      TTT4=SIGSTR(4)-SSS4                                               
      TTT5=SIGSTR(5)-SSS5                                               
      TTT6=SIGSTR(6)-SSS6                                               
      SSS1=SIGSTR(1)                                                    
      SSS2=SIGSTR(2)                                                    
      SSS3=SIGSTR(3)                                                    
      SSS4=SIGSTR(4)                                                    
      SSS5=SIGSTR(5)                                                    
      SSS6=SIGSTR(6)                                                    
      WRITE (6,*) 'TOTAL STR 3  1 = ',TTT1                              
      WRITE (6,*) '             2 = ',TTT2                              
      WRITE (6,*) '             3 = ',TTT3                              
      WRITE (6,*) '             4 = ',TTT4                              
      WRITE (6,*) '             5 = ',TTT5                              
      WRITE (6,*) '             6 = ',TTT6                              
      DO 1251 I=1,KNG                                                   
      ZVXC(I)=ZVP(I)                                                  
 1251 CONTINUE                                                          
C     CALL NON-LOCAL PART SUBROUTINE                                    
      CALL STRNL(ETOT1)                                           
C     NON-LOCAL PS AND EWALD PARTS                                      
      DO 1400 IS=1,6                                                    
      SIGSTR(IS)=SIGSTR(IS)+SIGNL(IS)+SIGEWA(IS)                        
 1400 CONTINUE                                                          
      WRITE (6,*) 'SIGEWA       1 = ',SIGEWA(1)                         
      WRITE (6,*) 'SIGEWA       2 = ',SIGEWA(2)                         
      WRITE (6,*) 'SIGEWA       3 = ',SIGEWA(3)                         
      WRITE (6,*) 'SIGEWA       4 = ',SIGEWA(4)                         
      WRITE (6,*) 'SIGEWA       5 = ',SIGEWA(5)                         
      WRITE (6,*) 'SIGEWA       6 = ',SIGEWA(6)                         
      WRITE (6,*) 'SIGNL        1 = ',SIGNL(1)                          
      WRITE (6,*) 'SIGNL        2 = ',SIGNL(2)                          
      WRITE (6,*) 'SIGNL        3 = ',SIGNL(3)                          
      WRITE (6,*) 'SIGNL        4 = ',SIGNL(4)                          
      WRITE (6,*) 'SIGNL        5 = ',SIGNL(5)                          
      WRITE (6,*) 'SIGNL        6 = ',SIGNL(6)                          
      TTT1=SIGSTR(1)-SSS1                                               
      TTT2=SIGSTR(2)-SSS2                                               
      TTT3=SIGSTR(3)-SSS3                                               
      TTT4=SIGSTR(4)-SSS4                                               
      TTT5=SIGSTR(5)-SSS5                                               
      TTT6=SIGSTR(6)-SSS6                                               
      SSS1=SIGSTR(1)                                                    
      SSS2=SIGSTR(2)                                                    
      SSS3=SIGSTR(3)                                                    
      SSS4=SIGSTR(4)                                                    
      SSS5=SIGSTR(5)                                                    
      SSS6=SIGSTR(6)                                                    
      WRITE (6,*) 'TOTAL STR 4  1 = ',TTT1                              
      WRITE (6,*) '             2 = ',TTT2                              
      WRITE (6,*) '             3 = ',TTT3                              
      WRITE (6,*) '             4 = ',TTT4                              
      WRITE (6,*) '             5 = ',TTT5                              
      WRITE (6,*) '             6 = ',TTT6                              
      SIGSTR(1)=SIGSTR(1)-ETOT1*TOTCH/UNIVOL                            
      SIGSTR(4)=SIGSTR(4)-ETOT1*TOTCH/UNIVOL                            
      SIGSTR(6)=SIGSTR(6)-ETOT1*TOTCH/UNIVOL                            
      TTT1=SIGSTR(1)-SSS1                                               
      TTT2=SIGSTR(2)-SSS2                                               
      TTT3=SIGSTR(3)-SSS3                                               
      TTT4=SIGSTR(4)-SSS4                                               
      TTT5=SIGSTR(5)-SSS5                                               
      TTT6=SIGSTR(6)-SSS6                                               
      SSS1=SIGSTR(1)                                                    
      SSS2=SIGSTR(2)                                                    
      SSS3=SIGSTR(3)                                                    
      SSS4=SIGSTR(4)                                                    
      SSS5=SIGSTR(5)                                                    
      SSS6=SIGSTR(6)                                                    
      WRITE (6,*) 'TOTAL STR 5  1 = ',TTT1                              
      WRITE (6,*) '             2 = ',TTT2                              
      WRITE (6,*) '             3 = ',TTT3                              
      WRITE (6,*) '             4 = ',TTT4                              
      WRITE (6,*) '             5 = ',TTT5                              
      WRITE (6,*) '             6 = ',TTT6                              
      WRITE (6,*) 'TOTAL SUMM   1 = ',SSS1                              
      WRITE (6,*) '             2 = ',SSS2                              
      WRITE (6,*) '             3 = ',SSS3                              
      WRITE (6,*) '             4 = ',SSS4                              
      WRITE (6,*) '             5 = ',SSS5                              
      WRITE (6,*) '             6 = ',SSS6                              
      IF (NBZTYP.LE.1) GO TO 9000
      DENOM = 1.0D0/FLOAT(KOPR)                                         
      SSS(1,1)=SIGSTR(1)                                                
      SSS(2,2)=SIGSTR(4)                                                
      SSS(3,3)=SIGSTR(6)                                                
      SSS(1,2)=SIGSTR(2)                                                
      SSS(1,3)=SIGSTR(3)                                                
      SSS(2,3)=SIGSTR(5)                                                
      SSS(2,1)=SSS(1,2)                                                 
      SSS(3,1)=SSS(1,3)                                                 
      SSS(3,2)=SSS(2,3)                                                 
      DO 1210 I=1,3                                                     
      DO 1220 J=1,3                                                     
      TTT(I,J)=0.0D0                                                    
 1220 CONTINUE                                                          
 1210 CONTINUE                                                          
!XOCL SPREAD DO /IQ
      DO 2400 IOP = 1,KOPR                                            
      DO 2200 I=1,3                                                     
      DO 2300 J=1,3                                                     
      SXX=SSS(1,1)                                                      
      SXY=SSS(1,2)                                                      
      SXZ=SSS(1,3)                                                      
      SYY=SSS(2,2)                                                      
      SYZ=SSS(2,3)                                                      
      SZZ=SSS(3,3)                                                      
      TXX=OO(1,I,IOP)*SSS(I,J)*OO(1,J,IOP)                              
      TXY=OO(1,I,IOP)*SSS(I,J)*OO(2,J,IOP)                              
      TXZ=OO(1,I,IOP)*SSS(I,J)*OO(3,J,IOP)                              
      TYY=OO(2,I,IOP)*SSS(I,J)*OO(2,J,IOP)                              
      TYZ=OO(2,I,IOP)*SSS(I,J)*OO(3,J,IOP)                              
      TZZ=OO(3,I,IOP)*SSS(I,J)*OO(3,J,IOP)                              
      TTT(1,1)=TTT(1,1) + OO(1,I,IOP)*SSS(I,J)*OO(1,J,IOP)              
      TTT(1,2)=TTT(1,2) + OO(1,I,IOP)*SSS(I,J)*OO(2,J,IOP)              
      TTT(1,3)=TTT(1,3) + OO(1,I,IOP)*SSS(I,J)*OO(3,J,IOP)              
      TTT(2,2)=TTT(2,2) + OO(2,I,IOP)*SSS(I,J)*OO(2,J,IOP)              
      TTT(2,3)=TTT(2,3) + OO(2,I,IOP)*SSS(I,J)*OO(3,J,IOP)              
      TTT(3,3)=TTT(3,3) + OO(3,I,IOP)*SSS(I,J)*OO(3,J,IOP)              
 2300 CONTINUE                                                          
 2200 CONTINUE                                                          
 2400 CONTINUE                                                        
!XOCL END SPREAD SUM(TTT)
      SIGSTR(1)=TTT(1,1)*DENOM                                          
      SIGSTR(2)=TTT(1,2)*DENOM                                          
      SIGSTR(3)=TTT(1,3)*DENOM                                          
      SIGSTR(4)=TTT(2,2)*DENOM                                          
      SIGSTR(5)=TTT(2,3)*DENOM                                          
      SIGSTR(6)=TTT(3,3)*DENOM                                          
      SSS1=SIGSTR(1)                                                    
      SSS2=SIGSTR(2)                                                    
      SSS3=SIGSTR(3)                                                    
      SSS4=SIGSTR(4)                                                    
      SSS5=SIGSTR(5)                                                    
      SSS6=SIGSTR(6)                                                    
      WRITE (6,*) 'TOTAL SUMM OP1 = ',SSS1                              
      WRITE (6,*) '             2 = ',SSS2                              
      WRITE (6,*) '             3 = ',SSS3                              
      WRITE (6,*) '             4 = ',SSS4                              
      WRITE (6,*) '             5 = ',SSS5                              
      WRITE (6,*) '             6 = ',SSS6                              
 9000 CONTINUE
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                 
      SUBROUTINE STRNL(ETOT1)                                     
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                 
      IMPLICIT REAL(A-H,O-Y)                                          
      IMPLICIT COMPLEX(Z)                                            
      INCLUDE 'PACVPP'                                                  
C     VPP-PARALLEL
!XOCL SUBPROCESSOR PS(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IP=(PS,INDEX=1:KNV3,PART=BAND)
      DIMENSION SSS(KNG1,KNV3,KTYP,10),SS2(KNG1,KNV3,KTYP,9)
     &         ,SS3(KNG1,KNV3,KTYP,4),ZZZ(KNG1,KEG,KNV3)
     &         ,ZFC(KEG,KNV3,KATM,10),ZFC2(KEG,KNV3,KATM,6)
     &         ,OCCUU(KEG,KNV3),RAA(KNG1,KNV3) 
!XOCL LOCAL SSS(:,/IP,:,:),SS2(:,/IP,:,:),SS3(:,/IP,:)
!XOCL LOCAL ZZZ(:,:,/IP),ZFC(:,/IP,:,:),ZFC2(:,/IP,:,:)
!XOCL LOCAL OCCUU(:,/IP),RAA(:,/IP)
      EQUIVALENCE (SNL,SSS),(SNL2,SS2),(SNL3,SS3),(ZAJ,ZZZ)
     &           ,(ZFBB,ZFC),(ZFBB2,ZFC2),(OCCUP,OCCUU)
     &           ,(RAK,RAA)
C     STRESS CALCULATION FOR NON-LOCAL PSEUDOPOTENTIAL PART             
      KSTR=KEG*KNV3*KATM*6                                              
      CCC =2.D0*RVOL/(2.D0*PAI)**3/FLOAT(KV3)                           
      RUNI=2.0D0*CCC                                                    
      DO 1003 IS=1,6                                                    
      SIGNL(IS)=0.0D0                                                   
 1003 CONTINUE                                                          
C     VPP-PARALLEL START 1
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 1901 IK=1,KV3                                                  
      DO 1002 IBAN=NBD1,NBD2                                            
      DO 1902 IA=1,KATM                                                 
      ZFC2(IBAN,IK,IA,1)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,2)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,3)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,4)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,5)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,6)=DCMPLX(0.0D0,0.0D0)                            
 1902 CONTINUE                                                          
 1002 CONTINUE                                                          
 1901 CONTINUE                                                          
!XOCL END SPREAD
C!XOCL END PARALLEL
C                                                                       
C     VPP-PARALLEL 2
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 1010 IK=1,KV3                                                  
      DO 1000 IA=1,KATM                                                 
      RWS =WS(KFTYPE(IA))*UNIVOL                                        
      RWP =WP(KFTYPE(IA))*UNIVOL                                        
      CWL(1)= RWS
      CWL(2)= RWP
      CWL(3)= RWP
      CWL(4)= RWP
      IF (NLSPD(KFTYPE(IA)).EQ.1) THEN
          LNUM  = 4
      ELSE
          LNUM  = 9
          RWD =WD(KFTYPE(IA))*UNIVOL
          CWL(5)= RWD
          CWL(6)= RWD
          CWL(7)= RWD
          CWL(8)= RWD
          CWL(9)= RWD
      END IF
C
      DO 1022 L=1,LNUM 
      DO 1020 IBAN=NBD1,NBD2                                            
      CW=0.5D0*RUNI*OCCUU(IBAN,IK)                                      
      STMP=CW*CWL(L)*DREAL(DCONJG(ZFC(IBAN,IK,IA,L))
     &              *ZFC(IBAN,IK,IA,L))
      SIGNL(1)=SIGNL(1)-STMP                                            
      SIGNL(4)=SIGNL(4)-STMP                                            
      SIGNL(6)=SIGNL(6)-STMP                                            
 1020 CONTINUE                                                          
 1022 CONTINUE
 1000 CONTINUE                                                          
 1010 CONTINUE                                                          
!XOCL END SPREAD SUM(SIGNL)
C!XOCL END PARALLEL
      SSS1=SIGNL(1)                                                     
      SSS2=SIGNL(2)                                                     
      SSS3=SIGNL(3)                                                     
      SSS4=SIGNL(4)                                                     
      SSS5=SIGNL(5)                                                     
      SSS6=SIGNL(6)                                                     
      WRITE (6,*) 'SIGNL     1  1 = ',SSS1                              
      WRITE (6,*) '             2 = ',SSS2                              
      WRITE (6,*) '             3 = ',SSS3                              
      WRITE (6,*) '             4 = ',SSS4                              
      WRITE (6,*) '             5 = ',SSS5                              
      WRITE (6,*) '             6 = ',SSS6                              
C                                                                       
      SK1=0.0D0                                                         
      SK2=0.0D0                                                         
      SK3=0.0D0                                                         
      SK4=0.0D0                                                         
      SK5=0.0D0                                                         
      SK6=0.0D0                                                         
C     VPP-PARALLEL START 3
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 1110 IK=1,KV3                                                  
      AKX = VX(IK)                                                      
      AKY = VY(IK)                                                      
      AKZ = VZ(IK)                                                      
      DO 1100 IA=1,KATM                                                 
      DO 1120 IBAN=NBD1,NBD2                                            
      CW=RUNI*OCCUU(IBAN,IK)                                            
      DO 1130 I=1,IBA(IK)                                               
      I1  = NBASE(I,IK)                                                 
      ZTMP=ZZZ(I,IBAN,IK)*DCONJG(ZFM2(I1,IA))               
     &    *SS2(I,IK,KFTYPE(IA),1)*RAA(I,IK)                            
      ZFC2(IBAN,IK,IA,1)=ZFC2(IBAN,IK,IA,1)                             
     &                  -ZTMP*(AKX+GX(I1))*(AKX+GX(I1))                 
      ZFC2(IBAN,IK,IA,2)=ZFC2(IBAN,IK,IA,2)                             
     &                  -ZTMP*(AKX+GX(I1))*(AKY+GY(I1))                 
      ZFC2(IBAN,IK,IA,3)=ZFC2(IBAN,IK,IA,3)                             
     &                  -ZTMP*(AKX+GX(I1))*(AKZ+GZ(I1))                 
      ZFC2(IBAN,IK,IA,4)=ZFC2(IBAN,IK,IA,4)                             
     &                  -ZTMP*(AKY+GY(I1))*(AKY+GY(I1))                 
      ZFC2(IBAN,IK,IA,5)=ZFC2(IBAN,IK,IA,5)                             
     &                  -ZTMP*(AKY+GY(I1))*(AKZ+GZ(I1))                 
      ZFC2(IBAN,IK,IA,6)=ZFC2(IBAN,IK,IA,6)                             
     &                  -ZTMP*(AKZ+GZ(I1))*(AKZ+GZ(I1))                 
 1130 CONTINUE                                                          
      SIGNL(1)=SIGNL(1)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,1)*DCONJG(ZFC(IBAN,IK,IA,1)) )            
      SIGNL(2)=SIGNL(2)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,2)*DCONJG(ZFC(IBAN,IK,IA,1)) )            
      SIGNL(3)=SIGNL(3)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,3)*DCONJG(ZFC(IBAN,IK,IA,1)) )            
      SIGNL(4)=SIGNL(4)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,4)*DCONJG(ZFC(IBAN,IK,IA,1)) )            
      SIGNL(5)=SIGNL(5)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,5)*DCONJG(ZFC(IBAN,IK,IA,1)) )            
      SIGNL(6)=SIGNL(6)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,6)*DCONJG(ZFC(IBAN,IK,IA,1)) )            
 1120 CONTINUE                                                          
 1100 CONTINUE                                                          
 1110 CONTINUE                                                          
!XOCL END SPREAD SUM(SIGNL)
C!XOCL END PARALLEL
      TTT1=SIGNL(1)-SSS1                                                
      TTT2=SIGNL(2)-SSS2                                                
      TTT3=SIGNL(3)-SSS3                                                
      TTT4=SIGNL(4)-SSS4                                                
      TTT5=SIGNL(5)-SSS5                                                
      TTT6=SIGNL(6)-SSS6                                                
      SSS1=SIGNL(1)                                                     
      SSS2=SIGNL(2)                                                     
      SSS3=SIGNL(3)                                                     
      SSS4=SIGNL(4)                                                     
      SSS5=SIGNL(5)                                                     
      SSS6=SIGNL(6)                                                     
      WRITE (6,*) 'SIGNL     2  1 = ',TTT1                              
      WRITE (6,*) '             2 = ',TTT2                             
      WRITE (6,*) '             3 = ',TTT3                              
      WRITE (6,*) '             4 = ',TTT4                              
      WRITE (6,*) '             5 = ',TTT5                              
      WRITE (6,*) '             6 = ',TTT6                              
C                  4                                                     
      DO 1300 L =2,9                                                    
C     VPP-PARALLEL START 4
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 1912 IK=1,KV3                                                  
      DO 1913 IA=1,KATM                                                 
      IF (NLSPD(KFTYPE(IA)).EQ.1 .AND. L.GT.4)
     &    GO TO 1913
      DO 1911 IBAN=NBD1,NBD2                                            
      ZFC2(IBAN,IK,IA,1)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,2)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,3)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,4)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,5)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,6)=DCMPLX(0.0D0,0.0D0)                            
 1911 CONTINUE                                                          
 1913 CONTINUE                                                          
 1912 CONTINUE                                                          
!XOCL END SPREAD
C!XOCL END PARALLEL
C
C     VPP-PARALLEL START 5
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 1320 IK=1,KV3                                                  
      AKX = VX(IK)                                                      
      AKY = VY(IK)                                                      
      AKZ = VZ(IK)                                                      
      DO 1310 IA=1,KATM                                                 
      IF (NLSPD(KFTYPE(IA)).EQ.1 .AND. L.GT.4)
     &    GO TO 1310
      DO 1330 IBAN=NBD1,NBD2                                            
      CW=RUNI*OCCUU(IBAN,IK)                                            
      DO 1340 I=1,IBA(IK)                                               
      I1  = NBASE(I,IK)                                                 
      ZTMP=ZZZ(I,IBAN,IK)*DCONJG(ZFM2(I1,IA))               
     &    *( SS2(I,IK,KFTYPE(IA),L) )*RAA(I,IK)                        
      ZFC2(IBAN,IK,IA,1)=ZFC2(IBAN,IK,IA,1)-ZTMP                        
     &                  *(AKX+GX(I1))*(AKX+GX(I1))                      
      ZFC2(IBAN,IK,IA,2)=ZFC2(IBAN,IK,IA,2)-ZTMP                        
     &                  *(AKX+GX(I1))*(AKY+GY(I1))                      
      ZFC2(IBAN,IK,IA,3)=ZFC2(IBAN,IK,IA,3)-ZTMP                        
     &                  *(AKX+GX(I1))*(AKZ+GZ(I1))                      
      ZFC2(IBAN,IK,IA,4)=ZFC2(IBAN,IK,IA,4)-ZTMP                        
     &                  *(AKY+GY(I1))*(AKY+GY(I1))                      
      ZFC2(IBAN,IK,IA,5)=ZFC2(IBAN,IK,IA,5)-ZTMP                        
     &                  *(AKY+GY(I1))*(AKZ+GZ(I1))                      
      ZFC2(IBAN,IK,IA,6)=ZFC2(IBAN,IK,IA,6)-ZTMP                        
     &                  *(AKZ+GZ(I1))*(AKZ+GZ(I1))                      
 1340 CONTINUE                                                          
      SIGNL(1)=SIGNL(1)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,1)*DCONJG( ZFC(IBAN,IK,IA,L)) )           
      SIGNL(2)=SIGNL(2)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,2)*DCONJG( ZFC(IBAN,IK,IA,L)) )           
      SIGNL(3)=SIGNL(3)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,3)*DCONJG( ZFC(IBAN,IK,IA,L)) )           
      SIGNL(4)=SIGNL(4)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,4)*DCONJG( ZFC(IBAN,IK,IA,L)) )           
      SIGNL(5)=SIGNL(5)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,5)*DCONJG( ZFC(IBAN,IK,IA,L)) )           
      SIGNL(6)=SIGNL(6)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,6)*DCONJG( ZFC(IBAN,IK,IA,L)) )           
 1330 CONTINUE                                                          
 1310 CONTINUE                                                          
 1320 CONTINUE                                                          
!XOCL END SPREAD SUM(SIGNL)
C!XOCL END PARALLEL
 1300 CONTINUE                                                          
      TTT1=SIGNL(1)-SSS1                                                
      TTT2=SIGNL(2)-SSS2                                                
      TTT3=SIGNL(3)-SSS3                                                
      TTT4=SIGNL(4)-SSS4                                                
      TTT5=SIGNL(5)-SSS5                                                
      TTT6=SIGNL(6)-SSS6                                                
      SSS1=SIGNL(1)                                                     
      SSS2=SIGNL(2)                                                     
      SSS3=SIGNL(3)                                                     
      SSS4=SIGNL(4)                                                     
      SSS5=SIGNL(5)                                                     
      SSS6=SIGNL(6)                                                     
      WRITE (6,*) 'SIGNL     3  1 = ',TTT1                              
      WRITE (6,*) '             2 = ',TTT2                              
      WRITE (6,*) '             3 = ',TTT3                              
      WRITE (6,*) '             4 = ',TTT4                              
      WRITE (6,*) '             5 = ',TTT5                              
      WRITE (6,*) '             6 = ',TTT6                              
      CWL(1)=1.0D0
      CWL(2)=1.0D0
      CWL(3)=1.0D0
      CWL(4)=1.0D0
      CWL(5)=2.0D0
      CWL(6)=2.0D0
      CWL(7)=2.0D0
      CWL(8)=2.0D0
      CWL(9)=2.0D0
      CWL(10)=1.0D0
C                 4                                                      
      DO 1400 L=2,10
C     VPP-PARALLEL START 6
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 1922 IK=1,KV3                                                  
      DO 1923 IA=1,KATM                                                 
      IF (NLSPD(KFTYPE(IA)).EQ.1 .AND. L.GT.4)
     &    GO TO 1923
      DO 1921 IBAN=NBD1,NBD2                                            
      ZFC2(IBAN,IK,IA,1)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,2)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,3)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,4)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,5)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,6)=DCMPLX(0.0D0,0.0D0)                            
 1921 CONTINUE                                                          
 1923 CONTINUE                                                          
 1922 CONTINUE                                                          
!XOCL END SPREAD
C!XOCL END PARALLEL
C
C     VPP-PARALLEL START 7
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 1420 IK=1,KV3                                                  
      AKX = VX(IK)                                                      
      AKY = VY(IK)                                                      
      AKZ = VZ(IK)                                                      
      DO 1410 IA=1,KATM                                                 
      IF (NLSPD(KFTYPE(IA)).EQ.1 .AND. L.GT.4)
     &    GO TO 1410
      DO 1430 IBAN=NBD1,NBD2                                            
      CW=CWL(L)*RUNI*OCCUU(IBAN,IK)                                            
      DO 1440 I=1,IBA(IK)                                               
      I1  = NBASE(I,IK)                                                 
      ZTMP=ZZZ(I,IBAN,IK)*DCONJG(ZFM2(I1,IA))               
     &    *SSS(I,IK,KFTYPE(IA),L)*RAA(I,IK)*RAA(I,IK)                   
      ZFC2(IBAN,IK,IA,1)=ZFC2(IBAN,IK,IA,1)+ZTMP                        
     &                  *(AKX+GX(I1))*(AKX+GX(I1))                      
      ZFC2(IBAN,IK,IA,2)=ZFC2(IBAN,IK,IA,2)+ZTMP                        
     &                  *(AKX+GX(I1))*(AKY+GY(I1))                      
      ZFC2(IBAN,IK,IA,3)=ZFC2(IBAN,IK,IA,3)+ZTMP                        
     &                  *(AKX+GX(I1))*(AKZ+GZ(I1))                      
      ZFC2(IBAN,IK,IA,4)=ZFC2(IBAN,IK,IA,4)+ZTMP                        
     &                  *(AKY+GY(I1))*(AKY+GY(I1))                      
      ZFC2(IBAN,IK,IA,5)=ZFC2(IBAN,IK,IA,5)+ZTMP                        
     &                  *(AKY+GY(I1))*(AKZ+GZ(I1))                      
      ZFC2(IBAN,IK,IA,6)=ZFC2(IBAN,IK,IA,6)+ZTMP                        
     &                  *(AKZ+GZ(I1))*(AKZ+GZ(I1))                      
 1440 CONTINUE                                                          
      SIGNL(1)=SIGNL(1)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,1)*DCONJG(ZFC(IBAN,IK,IA,L)) )            
      SIGNL(2)=SIGNL(2)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,2)*DCONJG(ZFC(IBAN,IK,IA,L)) )            
      SIGNL(3)=SIGNL(3)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,3)*DCONJG(ZFC(IBAN,IK,IA,L)) )            
      SIGNL(4)=SIGNL(4)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,4)*DCONJG(ZFC(IBAN,IK,IA,L)) )            
      SIGNL(5)=SIGNL(5)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,5)*DCONJG(ZFC(IBAN,IK,IA,L)) )            
      SIGNL(6)=SIGNL(6)+CW  *DREAL(                                     
     &        ZFC2(IBAN,IK,IA,6)*DCONJG(ZFC(IBAN,IK,IA,L)) )            
 1430 CONTINUE                                                          
 1410 CONTINUE                                                          
 1420 CONTINUE                                                          
!XOCL END SPREAD SUM(SIGNL)
C!XOCL END PARALLEL
 1400 CONTINUE                                                          
      TTT1=SIGNL(1)-SSS1                                                
      TTT2=SIGNL(2)-SSS2                                                
      TTT3=SIGNL(3)-SSS3                                                
      TTT4=SIGNL(4)-SSS4                                                
      TTT5=SIGNL(5)-SSS5                                                
      TTT6=SIGNL(6)-SSS6                                                
      SSS1=SIGNL(1)                                                     
      SSS2=SIGNL(2)                                                     
      SSS3=SIGNL(3)                                                     
      SSS4=SIGNL(4)                                                     
      SSS5=SIGNL(5)                                                     
      SSS6=SIGNL(6)                                                     
      WRITE (6,*) 'SIGNL     4  1 = ',TTT1                              
      WRITE (6,*) '             2 = ',TTT2                              
      WRITE (6,*) '             3 = ',TTT3                              
      WRITE (6,*) '             4 = ',TTT4                              
      WRITE (6,*) '             5 = ',TTT5                              
      WRITE (6,*) '             6 = ',TTT6                              
C                 1                                                      
      DO 2331 L=1,4
C     VPP-PARALLEL START 8
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 1942 IK=1,KV3                                                  
      DO 1943 IA=1,KATM                                                 
      IF (NLSPD(KFTYPE(IA)).EQ.1 .AND. L.GT.1)
     &    GO TO 1943
      DO 1941 IBAN=NBD1,NBD2                                            
      ZFC2(IBAN,IK,IA,1)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,2)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,3)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,4)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,5)=DCMPLX(0.0D0,0.0D0)                            
      ZFC2(IBAN,IK,IA,6)=DCMPLX(0.0D0,0.0D0)                            
 1941 CONTINUE                                                          
 1943 CONTINUE                                                          
 1942 CONTINUE                                                          
!XOCL END SPREAD
C!XOCL END PARALLEL
C
C     VPP-PARALLEL START 9
C!XOCL PARALLEL REGION
!XOCL SPREAD DO /IP
      DO 1520 IK=1,KV3                                                  
C     IWRT(IK) =IK                                                      
      AKX = VX(IK)                                                      
      AKY = VY(IK)                                                      
      AKZ = VZ(IK)                                                      
      DO 1510 IA=1,KATM                                                 
      IF (NLSPD(KFTYPE(IA)).EQ.1 .AND. L.GT.1)
     &    GO TO 1510
      IF (L.EQ.1) THEN
          RWC=1.0D0/(WP(KFTYPE(IA))*UNIVOL)                                
      ELSE
          RWC=1.0D0/(WD(KFTYPE(IA))*UNIVOL)                                
      END IF        
      DO 1530 IBAN=NBD1,NBD2                                            
      CW=RUNI*OCCUU(IBAN,IK)                                            
      DO 1540 I=1,IBA(IK)                                               
      I1  = NBASE(I,IK)                                                 
      ZTMP=DCONJG(ZZZ(I,IBAN,IK))*ZFM2(I1,IA)               
     &    *SS3(I,IK,KFTYPE(IA),L)*RAA(I,IK)                              
      ZFC2(IBAN,IK,IA,1)=ZFC2(IBAN,IK,IA,1)-ZTMP                        
     &                  *(AKX+GX(I1))                                   
      ZFC2(IBAN,IK,IA,2)=ZFC2(IBAN,IK,IA,2)-ZTMP                        
     &                  *(AKY+GY(I1))                                   
      ZFC2(IBAN,IK,IA,3)=ZFC2(IBAN,IK,IA,3)-ZTMP                        
     &                  *(AKZ+GZ(I1))                                   
 1540 CONTINUE                                                          
      SIGNL(1)=SIGNL(1)-RWC*CW  *DREAL(                                 
     &        ZFC2(IBAN,IK,IA,1)*DCONJG(ZFC2(IBAN,IK,IA,1)) )           
      SIGNL(2)=SIGNL(2)-RWC*CW  *DREAL(                                 
     &        ZFC2(IBAN,IK,IA,1)*DCONJG(ZFC2(IBAN,IK,IA,2)) )           
      SIGNL(3)=SIGNL(3)-RWC*CW  *DREAL(                                 
     &        ZFC2(IBAN,IK,IA,1)*DCONJG(ZFC2(IBAN,IK,IA,3)) )           
      SIGNL(4)=SIGNL(4)-RWC*CW  *DREAL(                                 
     &        ZFC2(IBAN,IK,IA,2)*DCONJG(ZFC2(IBAN,IK,IA,2)) )           
      SIGNL(5)=SIGNL(5)-RWC*CW  *DREAL(                                 
     &        ZFC2(IBAN,IK,IA,2)*DCONJG(ZFC2(IBAN,IK,IA,3)) )           
      SIGNL(6)=SIGNL(6)-RWC*CW  *DREAL(                                 
     &        ZFC2(IBAN,IK,IA,3)*DCONJG(ZFC2(IBAN,IK,IA,3)) )           
 1530 CONTINUE                                                          
 1510 CONTINUE                                                          
 1520 CONTINUE                                                          
!XOCL END SPREAD SUM(SIGNL)
C!XOCL END PARALLEL
 2331 CONTINUE
      TTT1=SIGNL(1)-SSS1                                                
      TTT2=SIGNL(2)-SSS2                                                
      TTT3=SIGNL(3)-SSS3                                                
      TTT4=SIGNL(4)-SSS4                                                
      TTT5=SIGNL(5)-SSS5                                                
      TTT6=SIGNL(6)-SSS6                                                
      SSS1=SIGNL(1)                                                     
      SSS2=SIGNL(2)                                                     
      SSS3=SIGNL(3)                                                     
      SSS4=SIGNL(4)                                                     
      SSS5=SIGNL(5)                                                     
      SSS6=SIGNL(6)                                                     
      WRITE (6,*) 'SIGNL     5  1 = ',TTT1                              
      WRITE (6,*) '             2 = ',TTT2                              
      WRITE (6,*) '             3 = ',TTT3                              
      WRITE (6,*) '             4 = ',TTT4                              
      WRITE (6,*) '             5 = ',TTT5                              
      WRITE (6,*) '             6 = ',TTT6                              
      RETURN                                                            
      END                                                               
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7  
      SUBROUTINE KSTEP                                                  
     >     (KNV3,NBZTYP,ALTV,RLTV,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,NFOUT,     
     <      KV3,VX,VY,VZ,QWGT)                                          
      IMPLICIT LOGICAL(A-Z)                                             
      INTEGER   KNV3,NBZTYP,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,KV3              
      INTEGER   INDX(0:20,0:20,0:20),IK,NFOUT                           
      REAL   ALTV(3,3),RLTV(3,3),A,B,C,PAI2,A1,B1,C1                  
      REAL   VX(KNV3),VY(KNV3),VZ(KNV3),QWGT(KNV3)                    
      LOGICAL  SPINP                                                    
      PAI2   = 8.D0*ATAN(1.D0)                                          
      SPINP  = .FALSE.                                                  
C***  IF((NBZTYP.EQ.0).OR.(NBZTYP.EQ.1)) THEN                           
      IF((NBZTYP.EQ.0).OR.(NBZTYP.EQ.1).OR.(NBZTYP.EQ.9)                
     > .OR.(NBZTYP.GE.11)) THEN
          IF(NBZTYP.EQ.0) THEN                                          
C----------------------------------------------                         
              WRITE(NFOUT,*) ' SURFACE B.Z. '                           
              CALL KPMSF                                                
     >         (KNV3,RLTV,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,NFOUT,             
     <          KV3,VX,VY,VZ)                                           
C----------------------------------------------                         
C****       ELSE IF(NBZTYP.EQ.1) THEN                                   
            ELSE IF(NBZTYP.EQ.1.OR.NBZTYP.EQ.9.OR.NBZTYP.GE.11) THEN
C----------------------------------------------                         
              WRITE(NFOUT,*) ' WHOLE B.Z. '                             
              CALL KPMWBZ                                               
     >         (KNV3,RLTV,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,NFOUT,             
     <          KV3,VX,VY,VZ)                                           
C----------------------------------------------                         
          END IF                                                        
          DO 100 IK=1,KNV3                                              
            QWGT(IK) = 1.D0/DFLOAT(KV3)                                 
  100     CONTINUE                                                      
        ELSE IF((NBZTYP.GE.2).AND.(NBZTYP.LE.5)) THEN                   
          IF(NBZTYP.EQ.2) THEN                                          
C----------------------------------------------------------             
              WRITE(NFOUT,*) ' SIMPLE CUBIC LATTICE '                   
              CALL SCCM                                                 
     >         (KNV3,NKX,SPINP,NFOUT,                                   
     <          KV3,VX,VY,VZ,QWGT,INDX)                                 
C----------------------------------------------------------             
            ELSE IF(NBZTYP.EQ.3) THEN                                   
C----------------------------------------------------------             
              WRITE(NFOUT,*) ' BCC LATTICE '                            
              CALL BCCM                                                 
     >         (KNV3,NKX,SPINP,NFOUT,                                   
     <          KV3,VX,VY,VZ,QWGT,INDX)                                 
C----------------------------------------------------------             
            ELSE IF((NBZTYP.EQ.4).OR.(NBZTYP.EQ.5)) THEN                
C----------------------------------------------------------             
              WRITE(NFOUT,*) ' FCC LATTICE '                            
              CALL FCCM                                                 
     >         (KNV3,NKX,SPINP,NFOUT,                                   
     <          KV3,VX,VY,VZ,QWGT,INDX)                                 
C----------------------------------------------------------             
          END IF                                                        
          A = ABS(ALTV(1,1)) * 2.D0                                     
          B      = PAI2/A                                               
          DO 200 IK = 1,KV3                                             
            VX(IK) = B*VX(IK)                                           
            VY(IK) = B*VY(IK)                                           
            VZ(IK) = B*VZ(IK)                                           
  200     CONTINUE                                                      
        ELSE IF(NBZTYP.EQ.6) THEN                                       
C----------------------------------------------------------             
          WRITE(NFOUT,*) ' HEX LATTICE '                                
          CALL  HEXM                                                    
     >     (KNV3,NKX,NKY,RLTV,SPINP,NFOUT,                              
     <      KV3,VX,VY,VZ,QWGT,INDX)                                     
C----------------------------------------------------------             
        ELSE IF (NBZTYP.EQ.8 .OR. NBZTYP.EQ.10) THEN                    
          WRITE (NFOUT,*) ' TETRAGONAL LATTICE '                        
          CALL TETRAH                                                   
     >     (KNV3,NKX,NKZ,RLTV,SPINP,NFOUT,                              
     <      KV3,VX,VY,VZ,QWGT,INDX)                                     
          A = ABS(ALTV(1,1))*2.0D0                                      
          A1     = PAI2/A                                               
          B = ABS(ALTV(3,3))*2.0D0                                      
          B1     = PAI2/B                                               
          WRITE (6,*) A,B,A1,B1,ALTV(1,1),ALTV(3,3)                     
          DO 201 IK = 1,KV3                                             
            VX(IK) = A1*VX(IK)                                          
            VY(IK) = A1*VY(IK)                                          
            VZ(IK) = B1*VZ(IK)                                          
  201     CONTINUE                                                      
        ELSE IF (NBZTYP.EQ.9.OR.NBZTYP.GE.11) THEN
          WRITE (NFOUT,*) ' ORTHORHONBIC LATTICE '                      
          CALL APBO2                                                    
     >     (KNV3,NKX,NKY,NKZ,RLTV,SPINP,NFOUT,                          
     <      KV3,VX,VY,VZ,QWGT,INDX)                                     
          A = ABS(ALTV(1,1))*2.0D0                                      
          A1     = PAI2/A                                               
          B = ABS(ALTV(2,2))*2.0D0                                      
          B1     = PAI2/B                                               
          C = ABS(ALTV(3,3))*2.0D0                                      
          C1     = PAI2/C                                               
          DO 202 IK = 1,KV3                                             
            VX(IK) = A1*VX(IK)                                          
            VY(IK) = B1*VY(IK)                                          
            VZ(IK) = C1*VZ(IK)                                          
  202     CONTINUE                                                      
        ELSE                                                            
          WRITE(NFOUT,*) ' NBZTYP ERR IN KSTEP NBUB = ',NBZTYP          
          STOP                                                          
      END IF                                                            
C         KV3=1                                                         
C         A = ABS(ALTV(1,1))                                            
C         A1     = PAI2/A                                               
C         B = ABS(ALTV(3,3))                                            
C         B1     = PAI2/B                                               
C         WRITE (6,*) A,B,A1,B1,ALTV(1,1),ALTV(3,3)                     
C         DO 202 IK = 1,KV3                                             
C           VX(IK) = A1*VX(IK)                                          
C           VY(IK) = A1*VY(IK)                                          
C           VZ(IK) = B1*VZ(IK)                                          
C 202     CONTINUE                                                      
      WRITE(NFOUT,300) (IK,VX(IK),VY(IK),VZ(IK),QWGT(IK),IK=1,KV3)      
  300 FORMAT((' ','IK = ',I3,3F10.6,3X,'QWGT = ',F10.6))                
      RETURN                                                            
      END                                                               
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7  
      SUBROUTINE KPMWBZ                                                 
     >     (KNV3,RLTV,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,NFOUT,                 
     <      KV3,VX,VY,VZ)                                               
C************************************************************           
      IMPLICIT LOGICAL(A-Z)                                             
      INTEGER IKP,KNV3,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,KV3,MM,I1,I2,I3,NFOUT 
      REAL   RLTV(3,3),VX(KNV3),VY(KNV3),VZ(KNV3),FVX,FVY,FVZ         
      REAL   DGVX,DGVY,DGVZ,FV1,FV2,FV3                               
      FVX=DFLOAT(NKX*2)                                                 
      FVY=DFLOAT(NKY*2)                                                 
      FVZ=DFLOAT(NKZ*2)                                                 
      IKP=0
      IF (IKP.EQ.0) THEN
        MM=0                                                            
C        DGVX=0.5D0*(RLTV(1,1)/FVX+RLTV(1,2)/FVY+RLTV(1,3)/FVZ)
C        DGVY=0.5D0*(RLTV(2,1)/FVX+RLTV(2,2)/FVY+RLTV(2,3)/FVZ)
C        DGVZ=0.5D0*(RLTV(3,1)/FVX+RLTV(3,2)/FVY+RLTV(3,3)/FVZ)
        DGVX=0.0D0
        DGVY=0.0D0
        DGVZ=0.0D0
        DO 100 I1=1,NKX2                                                
                      FV1=DFLOAT(I1-NKX-1)/FVX                          
        DO 100 I2=1,NKY2                                                
                      FV2=DFLOAT(I2-NKY-1)/FVY                          
        DO 100 I3=1,NKZ2                                                
                      FV3=DFLOAT(I3-NKZ-1)/FVZ                          
        MM=MM+1                                                         
        VX(MM)=RLTV(1,1)*FV1 + RLTV(1,2)*FV2 + RLTV(1,3)*FV3 +DGVX      
        VY(MM)=RLTV(2,1)*FV1 + RLTV(2,2)*FV2 + RLTV(2,3)*FV3 +DGVY      
        VZ(MM)=RLTV(3,1)*FV1 + RLTV(3,2)*FV2 + RLTV(3,3)*FV3 +DGVZ      
  100 CONTINUE                                                          
        WRITE(NFOUT,*) 'NUMBER OF GENERATED K-POINTS=',MM               
        IF(MM.NE.KNV3) WRITE(NFOUT,*) '**WARNING MM SHOULD BE KNV3**'   
        KV3=MM                                                          
      ELSE IF (IKP.EQ.1) THEN
        VX(1)=0.0D0
        VY(1)=0.0D0
        VZ(1)=0.0D0
        KV3=1
      END IF
      RETURN                                                            
      END                                                               
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7  
      SUBROUTINE KPMSF                                                  
     >     (KNV3,RLTV,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,NFOUT,                 
     <      KV3,VX,VY,VZ)                                               
      IMPLICIT LOGICAL(A-Z)                                             
      INTEGER   KNV3,NKX,NKY,NKZ,NKX2,NKY2,NKZ2,KV3,MM,I2,I3,NFOUT      
      REAL   RLTV(3,3),VX(KNV3),VY(KNV3),VZ(KNV3),FVY,FVZ             
      REAL   DGVY,DGVZ,FV2,FV3                                        
      FVY=DFLOAT(NKY*2)                                                 
      FVZ=DFLOAT(NKZ*2)                                                 
      MM=0                                                              
      DGVY=0.5D0*(RLTV(2,2)/FVY+RLTV(2,3)/FVZ)                          
      DGVZ=0.5D0*(RLTV(3,2)/FVY+RLTV(3,3)/FVZ)                          
      DO 100 I2=1,NKY2                                                  
             FV2=DFLOAT(I2-NKY-1)/FVY                                   
      DO 100 I3=1,NKZ2                                                  
             FV3=DFLOAT(I3-NKZ-1)/FVZ                                   
             MM=MM+1                                                    
             VX(MM) = 0.D0                                              
             VY(MM) = RLTV(2,2)*FV2 + RLTV(2,3)*FV3 + DGVY              
             VZ(MM) = RLTV(3,2)*FV2 + RLTV(3,3)*FV3 + DGVZ              
  100 CONTINUE                                                          
      WRITE(NFOUT,*) 'NUMBER OF GENERATED K-POINTS=',MM                 
      IF(MM.NE.KNV3) WRITE(NFOUT,*) '**WARN MM SHOULD BE KNV3**'        
      KV3=MM                                                            
      RETURN                                                            
      END                                                               
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7  
      SUBROUTINE DNV
     >     (N,K,XK,
     <      YB)
C**********************************************************             
      IMPLICIT REAL(A-H,O-Y)
      IMPLICIT COMPLEX(Z)
      INCLUDE 'PACVPP'
      DIMENSION XK(K),YB(K)
      IF(N.EQ.0) THEN                                                   
      DO 100 I=1,K
            XQ1   = XK(I)
            YB(I) = DSJBES(0,XQ1)
  100     CONTINUE
        ELSE IF(N.EQ.1) THEN                                            
        DO 110 I=1,K
            XQ1   = X(I)
            YB(I) = DSJBES(1,XQ1)
  110     CONTINUE
        ELSE IF(N.EQ.2) THEN                                            
        DO 120 I=1,K
            XQ1   = X(I)
            YB(I) = DSJBES(2,XQ1)
  120     CONTINUE
        ELSE IF(N.EQ.3) THEN                                            
        DO 130 I=1,K
            XQ1   = X(I)
            YB(I) = DSJBES(3,XQ1)
  130     CONTINUE
        ELSE IF(N.EQ.4) THEN                                            
        DO 140 I=1,K
            XQ1   = X(I)
            YB(I) = DSJBES(4,XQ1)
  140     CONTINUE
      END IF                                                            
      RETURN                                                            
      END                                                               
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      SUBROUTINE BUBBLE(N,IX,IY,IZ,BX,BY,BZ,BR) 
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      INTEGER   IX(1),IY(1),IZ(1)
      REAL    BX(1),BY(1),BZ(1),BR(1)
C
      DO 262 I=1,N-1
       DO 270 J=I+1,N
         IF(BR(J).LT.BR(I)) THEN
            CR   =BR(I)
            CX   =BX(I)
            CY   =BY(I)
            CZ   =BZ(I)
            LX   =IX(I)
            LY   =IY(I)
            LZ   =IZ(I)
            BR(I)=BR(J)
            BX(I)=BX(J)
            BY(I)=BY(J)
            BZ(I)=BZ(J)
            IX(I)=IX(J)
            IY(I)=IY(J)
            IZ(I)=IZ(J)
            BR(J)=CR
            BX(J)=CX
            BY(J)=CY
            BZ(J)=CZ
            IX(J)=LX
            IY(J)=LY
            IZ(J)=LZ
         END IF
  270 CONTINUE
  262 CONTINUE
      RETURN
      END
[先頭][Copyright © 1998-2013][Top]