C######################################################################
C#                                                                    #
C# MEAM.f                                                             #
C#  1.Total energy calculation by Modified Embedded Atom Method(MEAM) #
C#  2.Energetics in uniform dilation and compression                  #
C#  3.Energetics in one-direction-dilation and compression            #
C#  4.Force calculation by Modified Embedded Atom Method(MEAM)        #
C#  5.Molecular Dynamics (MD)                                         #
C#  6.Rose's energy calculation                                       #
C#  7.MEAM parameter file convert from Baskes's "meamf" to "meamPF.*" #
C#  8.Infinite body calculation....                                   #
C#  9.Argument of Embedding Function F(rho) in equilibrium states     #
C# 10. ....                                                           #
C#            based on M.I.Baskes,Phys.Rev.B,46,2727(1992)            #
C#                                                                    #
C#                           created          by   Kunio TAKAHASHI    #
C#                                            at   1997.9.4           #
C#                           finally modified by   Kunio TAKAHASHI    #
C#                                            at   1999.8.17          #
C#                                                                    #
C#     >>>  I and EAMers are strongly grateful for                    #
C#          prof.M.I.Baskes's kind and helpful suggestions.  <<<      #
C#                                                                    #
C#  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                          #
C#  % Copyright (C) Kunio TAKAHASHI, 1999- %                          #
C#  %        All rights reseved.           %                          #
C#  %      takahak@mep.titech.ac.jp        %                          #
C#  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                          #
C#                                                                    #
C#  - When you publish data calculated by this program or             #
C#    any program made by modification of this program,               #
C#    you have to write it in "acknowledgment".                       #
C#  - When you want to use this program,                              #
C#    it is better for you to collaborate with K.Takahashi.           #
C#    because program is always being upgraded.                       #
C#                                                                    #
C#                                                                    #
C######################################################################
C <In your FORTRAN comlilling>
C When you want to change the requiered memory in your computer, 
C please replace following expression in each program unit
C ( main, subroutines, and functions ).
C      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C where each parameter has following meaning
C Maximum Number of Atoms                      : mna
C Maximun Number of Character of Material Name : mncmn
C Maximum Number of the Kind of Atoms          : mnka
C Smallest REAL*8 Number                       : sr8n
C.........
C These parameters have constant value.
C In this program, such parameters are written in small capital.
C You have better not to use the same name as a variable.
C Variables are witten in large capital.
C
C
C
C <COMMON constant>
C  COMMON   /PI/PI
C        REAL*8 PI                :3.1415.....
C <COMMON variables>
C COMMON /NOATOM/NOATOM
C        INTEGER NOATOM           :number of the atoms
C COMMON /POSITN/X,Y,Z
C        REAL*8 X(mna)            :position X of i-th atom
C        REAL*8 Y(mna)            :         Y
C        REAL*8 Z(mna)            :         Z
C COMMON /MATERL/MATERL
C        CHARACTER*17 MATERL(mna) :name of i-th atom
C COMMON /PSYDAN/PSYDAN
C        INTEGER PSYDAN(mna)      :psyudo-Atomic Number of i-th atom
C            introduced in order to decrease shared memory and to use
C            Si1, Si2, etc... as well as  Si. The psyudo-Atomic Number
C            PSYDAN() is same as the ordering number of data block in
C            meamPF.* which is MEAM Parameter File.
C            ex.) in case of X,Y,Z,Mat
C                            0,0,0,Cu      PSYDAN(1)=5
C                            1,0,0,Si      PSYDAN(2)=1
C                            0,1,0,Ge      PSYDAN(3)=22
C                            0,0,1,Si      PSYDAN(4)=1
C                            1,1,0,Ca      PSYDAN(5)=35
C                            1,0,1,Si      PSYDAN(6)=1
C                            0,1,1,Ca      PSYDAN(7)=35
C                            1,1,1,Si      PSYDAN(8)=1
C                            ........, 
C                            PSYDAN(?)=5  means Cu
C                            PSYDAN(?)=1  means Si
C                            PSYDAN(?)=22 means Ge
C                            PSYDAN(?)=35 means Ca
C                            ........(numbering depend of meamPF.*).
C            If you modify the ordering of data block in the meamPF.*
C            above number will be changed.
C COMMON /NAMEOA/NAMEOA
C COMMON /REFSTR/REFSTR
C COMMON /ZI/ZI
C COMMON /ATMNUM/ATMNUM
C COMMON /ATMWEI/ATMWEI
C COMMON /ALPHA/ALPHA
C COMMON /BETA/BETA0,BETA1,BETA2,BETA3
C COMMON /ESUBLM/ESUBLM
C COMMON /AI/AI
C COMMON /SPDF/T0,T1,T2,T3
C COMMON /RHO0/RHO0
C COMMON /REQUIL/REQUIL
C  CHARACTER*17 NAMEOA(39) Name of Atom :element name
C  CHARACTER*17 REFSTR(39) Ref. lattice :lattice structure for refference
C  INTEGER      ZI(39)     Z_i          :number of 1st nabours in ref.struct.
C  INTEGER      ATMNUM(39) Atomic Number:atomic number
C  REAL*8       ATMWEI(39) Atomic weight:atomic weight
C  REAL*8       ALPHA(39)  alpha        :Rose energy parameter
C  REAL*8       BETA0(39)  beta0        :Atomic electron density parameter
C  REAL*8       BETA1(39)  beta1        :
C  REAL*8       BETA2(39)  beta2        :
C  REAL*8       BETA3(39)  beta3        :
C  REAL*8       ESUBLM(39) E_subli.     :Sublimation energy E_i^0 (eV)
C  REAL*8       AI(39)     A_i          :Ai Embedding energy coeff.
C  REAL*8       T0(39)     t0           :s p d f weighting coeef.
C  REAL*8       T1(39)     t1           :
C  REAL*8       T2(39)     t2           :
C  REAL*8       T3(39)     t3           :
C  REAL*8       RHO0(39)   rho_0        :nomalizing e- density for Rose energy
C  REAL*8       REQUIL(39) R_i^0        :equilibrium nearest^neighbour dist.(m)
C COMMON /SCREEN/TOSCRN
C        INTEGER TOSCRN                 : Type of screening function 
C      =0           : no screening
C      =1 (:default): MEAM92(M.I.Baskes,Phys.Rev.B,46,2727(1992))
C      =2           : ... not defined yet ...
C              ...
C COMMON /SCRFAC/SCRFAC
C         REAL*8  SCRFAC(mna,mna)     : Screening factor
C               defined in subroutine CALSCR
C               I-th atom sees the contribution of the electron density 
C               by J-th atom. However, The contribution is shaded by 
C               other atoms. Then, the electron density which I-th atom sees
C               is SCRFAC(I,J) x {the contiribution by J-th atom}.
C COMMON /RHOEXP/RHOEXP
C        INTEGER RHOEXP            : Expansion type of rho
C      =0 (:default): eq.9a  shown in M.I.Baskes,Phys.Rev.B,46,2727(1992)
C      =1           : eq.9   shown in M.I.Baskes,Phys.Rev.B,46,2727(1992)
C      =2           : eq.4c  shown in Baskes&Johnson, "MEAM potential for
C                                     HCP metal", Modelling Simul.Mater. 
C                                     Sci.Eng.,2(1994)147-163
C              ...
C COMMON /MEAMPF/MEAMPF
C         CHARACTER*17 MEAMPF      : Parameter File of MEAM
C   =meamPF     : default
C   =meamPF.HCP : Baskes&Johnson,Modelling Simul.Mater.Sci.Eng.,2(1994)147
C   =meamPF.BAS : converted from "meamf" by Baskes
C
C
C
C
C KT format data file
C "KT format data file" includes (i) comment lines, 
C (ii) a number which specifies the type of screening function, 
C (iii) xyz-position data and (iv) name of atom.
C ex.)
C 1:C four Si atoms
C 2: 0
C 3: +0.0000000E+000 +0.0000000E+000 +0.0000000E+000 Si
C 4: +2.3500000E-010 +0.0000000E+000 +0.0000000E+000 Si
C 5: +0.5000000E-010 +2.0351597E-010 +0.0000000E+000 Si
C 6: +0.5000000E-010 +6.7838657E-011 +1.9187670E-010 Si
C ...
C
C
C
C
C Data format for reduction calculation
C 
C 1:C <four Si atoms for reconstruction calculation>
C 2:C Starting with "C","%","c",and"#" mean being comment. 
C 3: 1               0               .\PF\meamPF.YG
C 4: +0.0000000E+000 +0.0000000E+000 +0.0000000E+000 Si -7 1
C 5: +2.3500000E-010 +0.0000000E+000 +0.0000000E+000 Si -1 1
C 6: +0.5000000E-010 +2.0351597E-010 +0.0000000E+000 Si  0 1
C 7: +0.5000000E-010 +6.7838657E-011 +1.9187670E-010 Si  0 1
C
C In first non-comment line,
C 3: 1               0               .\PF\meamPF.YG
C    ~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~
C    TOSCRN          RHOEXP          MEAMPF
C These three parameters is as follows
C COMMON /SCREEN/TOSCRN:INTEGER: type of screening
C COMMON /RHOEXP/RHOEXP:INTEGER: type of expression of rho expansion
C COMMON /MEAMPF/MEAMPF:CHAR*17: file name of parameter file with path
C COMMON /SEL/SEL(mna):INTEGER: couting factor of energy in ESELCT()
C                                SEL()=0 ; not counted
C                                SEL()=1 ; counted
C
C And next,
C 7: +0.5000000E-010 +6.7838657E-011 +1.9187670E-010 Si  0 1
C ...~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~ ~~ ~~ ~
C There are 5 paremeters ;
C X(mna)  Y(mna)  Z(mna)  MATERL(mna)  REFATM(mna) SEL(mna).
C X(),Y(),Z(), and MATERL() are mentioned above.
C REFATM() is introduced for reconstraction calculation.
C i.e.,
C COMMON /REFATM/REFATM
C INTEGER REFATM(mna)
C
C REFATM(i)= ..   | : positive number is of refering atom. 
C          =  3   |     i-th atom is moved corresponding REFATM(i)-th atom.
C          =  2   |     i.e., Relative position is conservative between 
C          =  1  -      i-th atom and REFATM(i)-th atom
C          =  0     : Fixed                           <-
C          = -1     : fixed in xy-plane    : moved in z-direction <-
C          = -2     : fixed in yz-plane    : moved in x-direction
C          = -3     : fixed in zx-plane    : moved in y-direction
C          = -4     : fixed in x-direction : moved in yz-plane
C          = -5     : fixed in y-direction : moved in zx-plane
C          = -6     : fixed in z-direction : moved in xy-plane    <-
C          = -7     : Free                           <-
C          = -8  -
C          = -9   | : not defined yet.
C            ...  |
C
C You may frequently use 0, -7, -1, and -6.
C
C
C
C
C <parameter file of MEAM calculation>
C A parameter file ( meamPF, meamPF.yg, etc...) is essential. 
C Full path or Relative path from program have to be set 
C in the variable MEAMPF.
C COMMON /MEAMPF/MEAMPF:CHAR*17 
C
C
C
C
C <Instruction>
C How to use subroutines in this program in order to make your own program.
C Subroutine TEMEAM is the basic procedure. 
C C#######################
C      SUBROUTINE TEMEAM
C C#######################
C       REAL*8 ETOT
C       INTEGER FLAG
C C --- set NOATOM,X,Y,Z,MATERL,TOSCRN---
C C     !!Caution!! TOSCRN may be changed in this subroutine 
C       CALL READKT                                       ......(1)
C C --- calculate Rij ---
C       CALL CALRIJ                                       ......(2)
C C --- calculate screening factor ---
C       CALL CALSCR                                       ......(3)
C C --- read meamPF.* ---
C       CALL SETATP(FLAG)                                 ......(4)
C       IF (FLAG.EQ.1) GOTO 999
C C --- calculate Etotal ---
C       CALL ETOTAL(ETOT)                                 ......(5)
C       WRITE(*,*) '  ################################'
C       WRITE(*,*) '  Etotal = ',ETOT,'(eV)'
C       WRITE(*,*) '  ################################'
C   999 RETURN
C       END
C Please see above subroutine. It is a sample. 
C The actual subroutine TEMEAM might be changed in the present version.
C (1) 
C    by this subroutine, following variables are set.
C    COMMON /SCREEN/TOSCRN     :INTEGER: type of screening function.
C    COMMON /POSITN/X(),Y(),Z():REAL*8 : xyz position of I-th atom
C    COMMON /MATERL/MATERL()   :CHAR*17: name of I-th atom
C    COMMON /NOATOM/NOATOM     :INTEGER: total number of the atoms
C (2) 
C    by this subroutine, following variables are set.
C    COMMON /RIJ/RIJ(i,j)        :REAL*8:distance 
C                                 between i-th atom and j-th atom
C    COMMON /XIJALF/XIJALF(i,j,k):REAL*8:cosine to x(k=1),y(k=2),z(k=3)
C                                 of vector from i-th atom to j-th atom
C (3) 
C    by this subroutine, following variables are set.
C    COMMON /SCRFAC/SCRFAC(i,j):REAL*8:Screening factor
C           I-th atom sees the contribution of the electron density 
C           by J-th atom. However, The contribution is shaded by 
C           other atoms. 
C           Then, the electron density which I-th atom sees is 
C           summation of SCRFAC(I,J) x {the contiribution by J-th atom}.
C    this subroutine uses TOSCRN:INTEGER:type of screening function 
C          =0 : no screening
C          =1 : MEAM92(M.I.Baskes,Phys.Rev.B,46,2727(1992))
C          =2 : ... not defined yet ...
C    even when TOSCRN=0, this subroutine have to be called,
C    because SCRFAC() is all set to be 1.
C (4) 
C    - Calling READPF(), data of parameter file are read in variables.
C    - Give psydo- atomic numbers to atoms read from data file.
C    by this subroutine, following variables are set.
C    COMMON /PSYDAN/PSYDAN(i):INTEGER:psyudo-Atomic Number of i-th atom
C (5) 
C    by this subroutine, 
C    ETOT:total energy is calculated as a summation of function EI(i),
C    where
C    E(i)=F(i)+ sigma PHI(i,j)
C                 j
C
C
C
C
C <1999.5.19>
C - Instruction was added in comment lines.
C - Speed-upped. 
C   PHI() and RHOIL() skipps summations when SCRFAC()=0.
C   This is an idea by Mr.C.Nara.
C
C
C
C
C######################################################################
      PROGRAM MEAM
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
      CHARACTER MENU*80
C
      REAL*8 PI
      COMMON /PI/PI
C
      INTEGER TOSCRN,RHOEXP
      CHARACTER*17 MEAMPF
      COMMON /SCREEN/TOSCRN /RHOEXP/RHOEXP
      COMMON /MEAMPF/MEAMPF
C
C --- initialize constant ---
      PI=3.14159265359D0
C --- default calculation type ---
      TOSCRN=1
      RHOEXP=0
      MEAMPF='meamPF'
C --- menu ---
  001 WRITE(*,*) 
      WRITE(*,*) '  ###################################'
      WRITE(*,*) '  Modified Embedded Atom Method      '
      WRITE(*,*) '                   by Kunio TAKAHASHI'
      WRITE(*,*) '  M.I.Baskes,Phys.Rev.B,46,2727(1992)'
      WRITE(*,*) '  ###################################'
      WRITE(*,*) '  STATUS OF CALCULTION'
      WRITE(*,'(2A,I3)')  '  Type of screening function : '
     & ,'default =  1 (MEAM92)       :',TOSCRN
      WRITE(*,'(2A,I3)')  '  Type of rho expansion      : '
     & ,'default =  0 (MEAM92)       :',RHOEXP
      WRITE(*,'(2A,A17)') '  Parameter File Name        : '
     & ,'default = "meamPF" (MEAM92) :',MEAMPF
  002 WRITE(*,*) '  MENU'
      WRITE(*,*) '  1 :Total energy calculation                       '
      WRITE(*,*) '  2 :Energetics in uniform dilation&compression     '
      WRITE(*,*) '  3 :Make scalable position data file for energetics'
      WRITE(*,*) '  4 :Rose energy calculation vs distance            '
      WRITE(*,*) '  5 :Argument of F(rho) in equilibrium states       '
      WRITE(*,*) '  6 :Distribution of rho as a fuction of distance   '
      WRITE(*,*) '  7 :Check of calculation process 2 (above)         '
      WRITE(*,*) '  8 :Energy v.s. x,y,z deviation of a atom          '
      WRITE(*,*) '  9 :Energy v.s. 1,2,3 deviation of a atom          '
      WRITE(*,*) ' 10 :Reconstruction calculation                     '
      WRITE(*,*) ' 11 :Param.convert from "meamf" to "meamPF.*"(plan) '
      WRITE(*,*) ' 12 :  ...not defined yet ...                       '
      WRITE(*,*) '  A :Change screening function                      '
      WRITE(*,*) '  B :Change expansion function of rho               '
      WRITE(*,*) '  C :Change MEAM Parameter File                     '
      WRITE(*,*) '  Q :Quit                                           '
      WRITE(*,*) '  Input 1,2,3,..., or Q'
C
C --- call selected operation ---
      READ (*,'(A)') MENU
      IF ((MENU(1:1).EQ.'#').OR.(MENU(1:1).EQ.'=')) THEN
        GOTO 001
      ELSEIF (MENU.EQ.'1')  THEN
        CALL TEMEAM
      ELSEIF (MENU.EQ.'2')  THEN
        CALL ENEGET
      ELSEIF (MENU.EQ.'3')  THEN
        CALL SPDMAK
      ELSEIF (MENU.EQ.'4')  THEN
        CALL ROSEEP
      ELSEIF (MENU.EQ.'5')  THEN
        CALL ARGOFF
      ELSEIF (MENU.EQ.'6')  THEN
        CALL RHOVSR
      ELSEIF (MENU.EQ.'7')  THEN
        CALL CALCHK
      ELSEIF (MENU.EQ.'8')  THEN
        CALL XYZDEV
      ELSEIF (MENU.EQ.'9')  THEN
        CALL ARBDEV
      ELSEIF (MENU.EQ.'10')  THEN
        CALL RECNST
      ELSEIF ((MENU.EQ.'A').OR.(MENU.EQ.'a')) THEN
        CALL SCRTYP
      ELSEIF ((MENU.EQ.'B').OR.(MENU.EQ.'b')) THEN
        CALL RHOTYP
      ELSEIF ((MENU.EQ.'C').OR.(MENU.EQ.'c')) THEN
        CALL SLCTPF
      ELSEIF ((MENU.EQ.'Q').OR.(MENU.EQ.'q')) THEN
        GOTO 999
      ELSE
        GOTO 002
      ENDIF
      WRITE(*,*) 
      WRITE(*,*) 
      WRITE(*,*) 
C     &'  Your operation is finished ? '
C     &'  Input any character to continue...'
C      READ (*,*) MENU
      GOTO 001
C
  999 STOP
      END
C
C
C
C
C
C
C
C
C
C DDDDD   OOO      N    N  OOO  TTTTTTT 
C D    D O   O     NN   N O   O    T                  !!!!!!!
C D    D O   O     N N  N O   O    T                  !     !
C D    D O   O     N  N N O   O    T                  !     !
C D    D O   O     N   NN O   O    T                  !    !
C DDDDD   OOO      N    N  OOO     T                 !     !
C                                                    !    !
C   M       M           D  I     F                   !    !
C   MM     MM           D       F                    !   !
C   M M   M M  OOO   DDDD  I  FFFFF Y     Y         !   !
C   M  M M  M O   O D   D  I    F    Y   Y          !   !
C   M   M   M O   O D   D  I    F     Y Y           !  !
C   M       M  OOO   DDDD  I    F      Y            !!!!
C                                     Y
C      B              L              Y
C      B              L                           !!!!!
C      BBBBB   EEEE   L   OOO  W     W     W     !     !
C      B    B EEEEEE  L  O   O  W   W W   W     !       !
C      B    B E       L  O   O   W W   W W       !     !
C      BBBBB   EEEE   L   OOO     W     W         !!!!!
C
C
C
C
C  < To: Students in Takahashi lab. >
C  Never modify following subroutines !
C  Please read above <Instruction>, then you can understand 
C  how to use subroutines in this program.
C  Takahashi will frequently improve following subroutines.
C  You can replace following subroutines by modified version.
C  Howver, if you use these subroutines with somoe modification,
C  you can not safely replace them with improved version. 
C  Takahashi strongly recommend NOT to modify following subroutines.
C  Good Luck !
C
C
C
C
C######################################################################
      SUBROUTINE RECNST
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C Calculation of reconstruction
C (1)->(4)->+-->(2)-->(3)-->Calculation-->+-->return
C           |                             |
C           +-<--<--<--<--<--<--<--<--<--<+
C           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~roop
C - read XYZMR data file
C -
C
      INTEGER TOSCRN,RHOEXP,NOATOM,FLAG
      CHARACTER*17 MEAMPF
      CHARACTER*17 MATERL(mna)
      COMMON /SCREEN/TOSCRN
      COMMON /RHOEXP/RHOEXP
      COMMON /MEAMPF/MEAMPF
      COMMON /NOATOM/NOATOM
      COMMON /MATERL/MATERL
C
C
C
C
      REAL*8 DSTEP,E,EMEMO,ECHK
      INTEGER ROOPC,EFLAG
      CHARACTER*32 FILINP,FILOUT,FOMODE
C
C --- set roop counter = 0 ---
      ROOPC=0
C --- announcement ---
      WRITE(*,*) 
      WRITE(*,*) 
      WRITE(*,*) ' ##############################'
      WRITE(*,*) ' #   Reconstruction of atom   #'
      WRITE(*,*) ' ##############################'
C     --- caution ---
      WRITE(*,*) ' !!!!!!! Do not forget !!!!!!!'
      WRITE(*,*) ' This command may change following parameters !'
      WRITE(*,*) '     TOSCRN (type of screening function)'
      WRITE(*,*) '     RHOEXP (type of rho expansion)'
      WRITE(*,*) '     MEAMPF (filename of MEAM parameter file)'
      WRITE(*,*) 
C
C --- read a XYZMR fomat file ---(1)
C     --- data file name ---
 0001 WRITE(*,*) ' File name to read'
      READ (*,'(A)') FILINP
C     --- check ---
      OPEN (1,FILE=FILINP,STATUS='OLD',ERR=0001)
      CLOSE(1)
C     --- read data file ---
      CALL RXYZMR(FILINP)
C
C --- create a XYZMR fomat file for output ---
C     --- data file name ---
 0011 WRITE(*,*) ' File name for output'
      READ (*,'(A)') FILOUT
C     --- check ---
      OPEN (1,FILE=FILOUT,STATUS='NEW',ERR=0011)
C --- select type of output ---
 0012 WRITE(*,*) ' output in Appending mode or OverWrite mode (A/O) ?'
      READ (*,'(A)') FOMODE
      IF (     (FOMODE.NE.'A').AND.(FOMODE.NE.'a')
     &    .AND.(FOMODE.NE.'O').AND.(FOMODE.NE.'o')) GOTO 0012
      IF ((FOMODE(1:1).EQ.'A').OR.(FOMODE(1:1).EQ.'a')) THEN
        CALL OUTHD1(1,FILINP)
      ENDIF
      CLOSE(1)
C
C --- read meamPF.* ---(4)
      CALL SETATP(FLAG)
      IF (FLAG.EQ.1) THEN
        WRITE(*,*) ' failed in reading parameter file (SETATP)'
        GOTO 9999
      ENDIF
C
C --- input DSTEP ---
 0003 WRITE(*,*) 
     &' Input step of differentials',
     &' (1.e-13 to 1.e-20, larger maybe conbenient...)'
      READ (*,*,ERR=0003) DSTEP
C
C --- calculate initial energy (for output) ---
C     --- calculate Rij ---(2)
C      WRITE(*,*) ' calling CALRIJ...'
      CALL CALRIJ
C     --- calculate screening factor ---(3)
C      WRITE(*,*) ' calling CALSCR...'
      CALL CALSCR
C     --- calculate initial energy ---
C      WRITE(*,*) ' calling ESELCT...'
      CALL ESELCT(E)
      WRITE(*,*) ' Initial Total Energy (eV) =',E
C
C --- roop start ---
      EFLAG=0
C     ======= write data start =======
 1000 OPEN (1,FILE=FILOUT,STATUS='OLD')
C      WRITE(*,*) ' Data file opened to add...'
C     --- check ---
      IF ((FOMODE(1:1).EQ.'A').OR.(FOMODE(1:1).EQ.'a')) THEN
C       --- goto end of file ---
        CALL GO2EOF(1)
      ELSE
        CALL OUTHD1(1,FILINP)
      ENDIF
C     --- write header 2 ---
      CALL OUTHD2(1,ROOPC,E,DSTEP)
C     --- write x,y,z,m,r ---
      CALL WXYZMR(1)
C     --- close ---
      CLOSE(1)
C     ======= write data end =======
C     --- increase counter ---
      ROOPC=ROOPC+1
C     --- optimize on one direction in super space ---
      EMEMO=E
      CALL GO1WAY(DSTEP,EFLAG,E)
C     --- check quit or continue ---
      ECHK=DABS(E-EMEMO)
      IF (ECHK.LT.1.0D-250) THEN
        WRITE(*,*) ' cannot calculate precise energy more than now...'
        GOTO 9000
      ENDIF
      IF (EFLAG.EQ.1) THEN
        WRITE(*,*) ' eough precision in movement you requested '
        GOTO 9000
      ENDIF
C     ---  ---
      GOTO 1000
C --- roop end ---
C --- final output ---
 9000 OPEN (1,FILE=FILOUT,STATUS='OLD')
C      WRITE(*,*) ' Data file opened to add...'
C     --- check ---
      IF ((FOMODE(1:1).EQ.'A').OR.(FOMODE(1:1).EQ.'a')) THEN
C       --- goto end of file ---
        CALL GO2EOF(1)
      ELSE
        CALL OUTHD1(1,FILINP)
      ENDIF
C     --- write header 2 ---
      CALL OUTHD2(1,ROOPC,E,DSTEP)
C     --- write x,y,z,m,r ---
      CALL WXYZMR(1)
C     --- close ---
      CLOSE(1)
C
 9999 RETURN
      END
C################## subroutine only for RECNST ########################
      SUBROUTINE OUTHD2(N,R,E,D)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
      REAL*8 E,D
      INTEGER N,R
C
      WRITE(N,'(A)') '# '
      WRITE(N,'(A,I7,A)') 
     & '####################',R,'-th ###################'
      WRITE(N,'(A,1P,E25.16E3)') 
     & '# Etotal(of selected atoms)in eV =',E
      WRITE(N,'(A,1P,E25.16E3)') 
     & '# DSTEP in calculation           =',D
C
      RETURN
      END
C
C################## subroutine only for RECNST ########################
      SUBROUTINE OUTHD1(N,FILINP)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
      INTEGER TOSCRN,RHOEXP,NOATOM
      CHARACTER*17 MEAMPF
      COMMON /SCREEN/TOSCRN
      COMMON /RHOEXP/RHOEXP
      COMMON /MEAMPF/MEAMPF
      COMMON /NOATOM/NOATOM
C
      INTEGER N
      CHARACTER*32 FILINP
C
      WRITE(N,'(A)') '# '
      WRITE(N,'(A)') '# '
      WRITE(N,'(A)') '# Reconstruction Calculation by MEAM.f '
      WRITE(N,'(A)') '# '
      WRITE(N,'(A,I7)') 
     &'#   Type of screening function     (s)  :',TOSCRN
      WRITE(N,'(A,I7)') 
     &'#   Type of rho expansion          (r)  :',RHOEXP
      WRITE(N,'(A,A)') 
     &'#   MEAM parameter file            (p)  :',MEAMPF
      WRITE(N,'(A,A)') 
     &'#   Initial Data (s,r,p/x,y,z,m,r) file :',FILINP
      WRITE(N,'(A,I7)') 
     &'#   Number of atom                      :',NOATOM
      WRITE(N,'(A)') '# '
C
 9999 RETURN
      END
C######################################################################
      SUBROUTINE GO1WAY(DSTEP,EFLAG,E3)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C
      REAL*8 DSTEP
      REAL*8 E0,DIREX(mna),DIREY(mna),DIREZ(mna)
      REAL*8 EL,EC,ER,TL,TC,TR,TSTEP
      REAL*8 E1,E2,E3,E4,E5
      REAL*8 T1,T2,T3,T4,T5,DT2
      INTEGER EFLAG
      INTEGER I
      REAL*8 ESLONT
C      REAL*8 ECHK,EREF
C
      INTEGER NOATOM
      COMMON /NOATOM/NOATOM
      INTEGER NOMATM
      COMMON /NOMATM/NOMATM
      REAL*8 X(mna),Y(mna),Z(mna)
      COMMON /POSITN/X,Y,Z
C
C
C --- initiallize flag ---
      EFLAG=0
C --- calculate Rij ---
      CALL CALRIJ
C --- calculate screening factor ---
      CALL CALSCR
C --- calculate initial energy ---
      CALL ESELCT(E0)
C --- calculate direction (unit vector) ---
      WRITE(*,*) '  call DIRECT in GO1WAY: direction calculating...'
      CALL DIRECT(E0,DSTEP,DIREX,DIREY,DIREZ)
C
C
C --- if you want to accerelate calculation, select larger TSTEP. ---
C     if you want to avoid local minimum trap, select larger TSTEP.
C     TSTEP=DSTEP
      TSTEP=DSTEP*10
C     TSTEP=DSTEP*100
C     TSTEP=DSTEP*1000
C     TSTEP=DSTEP*10000
C --- although smaller TSETP is more safety... ---
C
C === go downward to find bottom with acceralation start ===
C     TL=-TSTEP,EL=+++(TL)
C     TC=0.D.  ,EC=E0
C     TR=TSTEP ,ER=+++(TR)
C     if E1>E0, very interesting ! multi body effect occores !
C        then, output data to file.
C     if EL<R0>ER, very interesting ! saddle point found !
C        then, output data to file.
C     --- roop1 ----
C     +--> if EL>EC<ER, find bottom -->-->-->-->--+
C     |    if EL>ER,                              |
C     |           TC<-TR(EC<-ER)                  |
C     |           TR=TC+(TC-TL)                   |
C     |           ER=+++(TR)                      |
C     |    if EL<ER,                              |
C     |           TC<-TL(EC<-EL)                  |
C     |           TL=TC-(TR-TC)                   |
C     |           EL=+++(TL)                      |
C     +-----------+                               |
C === go down with accerelation end ===  <--<--<--+
C
C === find bottom with breaking start ===
C     T1=TL,E1=EL
C     T3=TC,E3=EC
C     T5=TR,E5=ER
C     --- roop2 ----
C     +--> set T2, calculate E2=+++(T2)
C     |    set T4, calculate E4=+++(T4)
C     |    if E1>E2<E3 (E2 is minimum),
C     |           T5=T3,E5=E3
C     |           T3=T2,E3=E2
C     |    if E2>E3<E4 (E3 is minimum),
C     |           T1=T2,E1=E2
C     |           T5=T4,E5=E4
C     |    if E3>E4<E5 (E4 is minimum),
C     |           T1=T3,E1=E3
C     |           T3=T4,E3=E4
C     |    if (T4-T2)<DSTEP----------------+
C     +-----------+                        |
C === find bottom with breaking end === <--+
C
C === go downward to find bottom with acceralation ===
      WRITE(*,*) '  go downward to find bottom with acceralation...'
      TR=TSTEP
      ER=ESLONT(TR,DIREX,DIREY,DIREZ)
      TC=0
      EC=E0
      TL=-TSTEP
      EL=ESLONT(TL,DIREX,DIREY,DIREZ)
C     ---- you should detect following 2 interesting phenomena ---
C     if E1>E0, very interesting ! multi body effect occores !
C        then, output data to file.
C     if EL<R0>ER, very interesting ! saddle point found !
C        then, output data to file.
C     ---- however now pending ---
C     --- roop1 ----
  100 CONTINUE
C      WRITE(*,*) 'Tl,Tc,Tr',TL,TC,TR
C      WRITE(*,*) 'El,Ec,Er',EL,EC,ER
      IF ((EL.GT.EC).AND.(EC.LT.ER)) GOTO 199
      IF    ((EL.LT.EC).AND.(EC.GT.ER)) THEN
        WRITE(*,*) '  You found saddle point ! go to steep way..'
      ELSEIF (EL.GT.ER) THEN
        TC=TR
        EC=ER
        TR=TC*2-TL
        ER=ESLONT(TR,DIREX,DIREY,DIREZ)
      ELSEIF (EL.LT.ER) THEN
        TC=TL
        EC=EL
        TL=TC*2-TR
        EL=ESLONT(TL,DIREX,DIREY,DIREZ)
      ELSEIF ((EL.EQ.EC).AND.(EC.EQ.ER)) THEN
        WRITE(*,*) '  direction lost... strange... '
        GOTO 199
      ENDIF
      GOTO 100
  199 CONTINUE
C      WRITE(*,*) '   roop1 end in GO1WAY...'
C
C
C === find bottom with breaking start ===
      WRITE(*,*) '  finding bottom with breaking...'
      T1=TL
      E1=EL
      T3=TC
      E3=EC
      T5=TR
      E5=ER
C     --- roop2 ----
  200 CONTINUE
C      WRITE(*,*) 'T1,T3,T5',T1,T3,T5
C      WRITE(*,*) 'E1,E3,E5',E1,E3,E5
      T2=(T1+T3)/2
      E2=ESLONT(T2,DIREX,DIREY,DIREZ)
      T4=(T3+T5)/2
      E4=ESLONT(T4,DIREX,DIREY,DIREZ)
C     IF    ((E2.LT.E3).AND.(E2.LT.E4)) THEN
      IF    ((E2.LE.E3).AND.(E2.LE.E4)) THEN
        T5=T3
        E5=E3
        T3=T2
        E3=E2
C     ELSEIF((E3.LT.E2).AND.(E3.LT.E4)) THEN
      ELSEIF((E3.LE.E2).AND.(E3.LE.E4)) THEN
        T1=T2
        E1=E2
        T5=T4
        E5=E4
C     ELSEIF((E4.LT.E2).AND.(E4.LT.E3)) THEN
      ELSEIF((E4.LE.E2).AND.(E4.LE.E3)) THEN
        T1=T3
        E1=E3
        T3=T4
        E3=E4
      ENDIF
      DT2=T4-T2
      IF (DT2.GE.DSTEP) GOTO 200
C      WRITE(*,*) '   roop2 end in GO1WAY...'
C
C === move to T3 poinot in (DIREY(),DIREY(),DIREZ()) direction ---
      DO 002 I=1,NOATOM
        X(I)=X(I)+T3*DIREX(I)
        Y(I)=Y(I)+T3*DIREY(I)
        Z(I)=Z(I)+T3*DIREZ(I)
  002 CONTINUE
C
C      WRITE(*,'(A7,1P,E25.16E3,A26)') 
      WRITE(*,*) 
     & ' Moved ',T3,' (m) to (grad E) direction. '
C      WRITE(*,'(A23,1P,E25.16E3,A4)') 
      WRITE(*,*) 
     & ' Turning here ! Etotal= ',E3,' (eV) '
C
C
C === end flag check ===
C    --- differencial has meaning ? ---
      IF (DABS(T3).LE.DSTEP) EFLAG=1
CC    --- energy per atom has meaning ? ---
C     ECHK=DABS(E5-E3)+DABS(E1-E3)
C     EREF=1.D-5*E3/NOMATM
C     IF (ECHK.LT.0.E-3) EFLAG=2
C
C
      RETURN
      END
C======================================================
      FUNCTION ESLONT(T,EX,EY,EZ)
C======================================================
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C Calculation of total energy ( of selected atoms ) on t-axis
C only call from subroutine GO1WAY()
C EY(),EY(),EZ():unit vector of t-axis in super space 
C T             :displace ment from X()Y()Z().
C X()Y()Z() is conservative !
C 
      REAL*8 ESLONT,T,EX(mna),EY(mna),EZ(mna)
      REAL*8 XMEMO(mna),YMEMO(mna),ZMEMO(mna)
      COMMON /POSITN/X,Y,Z
      REAL*8 X(mna),Y(mna),Z(mna)
      INTEGER NOATOM
      COMMON /NOATOM/NOATOM
C
      REAL*8 E
      INTEGER I
C
C      WRITE(*,*) '  ESLONT called '
C --- back up positon data X()Y()Z() ---
      DO 001 I=1,NOATOM
        XMEMO(I)=X(I)
        YMEMO(I)=Y(I)
        ZMEMO(I)=Z(I)
  001 CONTINUE
C
C --- move T to (EY(),EY(),EZ()) direction ---
      DO 002 I=1,NOATOM
        X(I)=XMEMO(I)+T*EX(I)
        Y(I)=YMEMO(I)+T*EY(I)
        Z(I)=ZMEMO(I)+T*EZ(I)
  002 CONTINUE
C
C --- calculate Rij ---
      CALL CALRIJ
C --- calculate screening factor ---
      CALL CALSCR
C --- calculate initial energy ---
      CALL ESELCT(E)
      ESLONT=E
C
C --- recover positon data X()Y()Z() ---
      DO 003 I=1,NOATOM
        X(I)=XMEMO(I)
        Y(I)=YMEMO(I)
        Z(I)=ZMEMO(I)
  003 CONTINUE
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE DIRECT(E0,DSTEP,DIREX,DIREY,DIREZ)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C calculate direction to move
C giving
C   E0     :toal energy (selected) for given status
C   DSTEP  :differential step for numerical calculation
C this subrouting returns an unit vector to move : 
C   DIREX(),DIREY(),DIREZ()
C
      REAL*8 DSTEP,DIREX(mna),DIREY(mna),DIREZ(mna)
      REAL*8 E0,VLEN
      INTEGER I
C
      INTEGER NOATOM
      COMMON /NOATOM/NOATOM
      INTEGER REFATM(mna)
      COMMON /REFATM/REFATM
C
C --- initialize vector length ---
      VLEN=0.D0
C --- calculate (grad E) for REFATM()=<0 atoms ---
C      WRITE(*,*) '  calculating (grad E) for REFATM()=<0 atoms ... '
      DO 001 I=1,NOATOM
C        WRITE(*,*) '  i = ',I
        IF     (REFATM(I).EQ.0)  THEN
          DIREX(I)=0
          DIREY(I)=0
          DIREZ(I)=0
        ELSEIF (REFATM(I).EQ.-1) THEN
          CALL GRADE(E0,I,DSTEP,DIREX(I),DIREY(I),DIREZ(I))
          DIREX(I)=0
          DIREY(I)=0
          VLEN=VLEN+DIREZ(I)**2
        ELSEIF (REFATM(I).EQ.-2) THEN
          CALL GRADE(E0,I,DSTEP,DIREX(I),DIREY(I),DIREZ(I))
          DIREY(I)=0
          DIREZ(I)=0
          VLEN=VLEN+DIREX(I)**2
        ELSEIF (REFATM(I).EQ.-3) THEN
          CALL GRADE(E0,I,DSTEP,DIREX(I),DIREY(I),DIREZ(I))
          DIREZ(I)=0
          DIREX(I)=0
          VLEN=VLEN+DIREY(I)**2
        ELSEIF (REFATM(I).EQ.-4) THEN
          CALL GRADE(E0,I,DSTEP,DIREX(I),DIREY(I),DIREZ(I))
          DIREX(I)=0
          VLEN=VLEN+DIREY(I)**2+DIREZ(I)**2
        ELSEIF (REFATM(I).EQ.-5) THEN
          CALL GRADE(E0,I,DSTEP,DIREX(I),DIREY(I),DIREZ(I))
          DIREY(I)=0
          VLEN=VLEN+DIREX(I)**2+DIREZ(I)**2
        ELSEIF (REFATM(I).EQ.-6) THEN
          CALL GRADE(E0,I,DSTEP,DIREX(I),DIREY(I),DIREZ(I))
          DIREZ(I)=0
          VLEN=VLEN+DIREX(I)**2+DIREY(I)**2
        ELSEIF (REFATM(I).EQ.-7) THEN
          CALL GRADE(E0,I,DSTEP,DIREX(I),DIREY(I),DIREZ(I))
          VLEN=VLEN+DIREX(I)**2+DIREY(I)**2+DIREZ(I)**2
        ENDIF
C      WRITE(*,*) ' I,dXdYdZ',I,DIREX(I),DIREY(I),DIREZ(I)
  001 CONTINUE
C --- normalize (grad E) for REFATM()=<0 atoms ---
C      WRITE(*,*) '  normalizing (grad E) for REFATM()=<0 atoms ... '
      VLEN=DSQRT(VLEN)
      IF (VLEN.LT.1.D0) THEN
        WRITE(*,*) '    (grad Etot)<1(eV/m). It seems enough small !'
      ENDIF
      DO 002 I=1,NOATOM
        IF (REFATM(I).LE.0)  THEN
          DIREX(I)=DIREX(I)/VLEN
          DIREY(I)=DIREY(I)/VLEN
          DIREZ(I)=DIREZ(I)/VLEN
        ENDIF
  002 CONTINUE
C
C --- set (grad E) for refering (REFATM()>0) atoms ---
C      WRITE(*,*) '  setting (grad E) for refering (REFATM()>0) atoms...'
      DO 003 I=1,NOATOM
        IF     (REFATM(I).GT.NOATOM)  THEN
          WRITE(*,*) '  !!! WARNING !!!'
          WRITE(*,*) '  restriction not defined in DIRECT()'
          STOP
        ELSEIF (REFATM(I).GT.0) THEN
          DIREX(I)=DIREX(REFATM(I))
          DIREY(I)=DIREY(REFATM(I))
          DIREZ(I)=DIREZ(REFATM(I))
        ENDIF
  003 CONTINUE
C
      RETURN
      END
C######################################################################
      SUBROUTINE GRADE(E0,I,DSTEP,GX,GY,GZ)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C After steps (1), (2), (3), and (4), (see above instruction), 
C set E0    :REAL*8  :total energy of selected atoms
C     I     :INTEGER
C     DSTEP :REAL*8  :step for differential.
C and call this subroutine.
C Then it returns GX=dEtotal/dxi
C                 GY=dEtotal/dyi
C                 GZ=dEtotal/dzi
C by moving I-th atom.
C The position of atoms is consevative.
C Differential step DSTEP should be 1/10 of precision.
C
C
      REAL*8 DSTEP,GX,GY,GZ
      INTEGER I
      REAL*8 X(mna),Y(mna),Z(mna)
      COMMON /POSITN/X,Y,Z
C
      REAL*8 XMEMO,YMEMO,ZMEMO,E0,E1
C
C --- save the original position data ---
      XMEMO=X(I)
      YMEMO=Y(I)
      ZMEMO=Z(I)
CC --- calculate total energy for original position ---
C      this line have become not-needed by giving E0 as an argument
C      for the purpose of speed up
C      CALL ESELCT(E0)
C      however,
C      as shown following comment,
C     ---- energy associated with I-th atom should be calculated ----
C      ESELCT(E0) is replaced to anather subroutine
C      that subroutine should be called here...
C
C --- calculate GX ---
C      WRITE(*,*)  ' dEtotal/dxi'
C     --- deviate X(I) -> X(I)+DSTEP ---
      X(I)=XMEMO+DSTEP
C     ====<following step (2) and (3) should be modified>====
C     ---- parameters concerning I-th atom should be modified ----
C     --- calculate Rij step(2) ---
      CALL CALRIJ
C      CALL MODRIJ(I) : here CALRIJ must be called
C     --- calculate screening factor step(3) ---
      CALL CALSCR
C     ====<above step (2) and (3) should be modified>====
C     --- calculate total energy ---
C     ---- energy associated with I-th atom should be calculated ----
      CALL ESELCT(E1)
      GX=(E1-E0)/DSTEP
C
C --- calculate GY ---
C      WRITE(*,*)  ' dEtotal/dyi'
C     --- deviate Y(I) -> Y(I)+DSTEP ---
      Y(I)=YMEMO+DSTEP
C     ====<following step (2) and (3) should be modified>====
C     ---- parameters concerning I-th atom should be modified ----
C     --- calculate Rij step(2) ---
C      CALL CALRIJ
C      WRITE(*,*)  ' calling MODRIJ for dEtotal/dyi'
      CALL MODRIJ(I)
C     --- calculate screening factor step(3) ---
C      WRITE(*,*)  ' calling CALSCR for dEtotal/dyi'
      CALL CALSCR
C     ====<above step (2) and (3) should be modified>====
C     --- calculate total energy ---
C     ---- energy associated with I-th atom should be calculated ----
C      WRITE(*,*)  ' calling Eselect for dEtotal/dyi'
      CALL ESELCT(E1)
      GY=(E1-E0)/DSTEP
C
C --- calculate GZ ---
C      WRITE(*,*)  ' dEtotal/dzi'
C     --- deviate Z(I) -> Z(I)+DSTEP ---
      Z(I)=ZMEMO+DSTEP
C     ====<following step (2) and (3) should be modified>====
C     ---- parameters concerning I-th atom should be modified ----
C     --- calculate Rij step(2) ---
C      CALL CALRIJ
      CALL MODRIJ(I)
C     --- calculate screening factor step(3) ---
      CALL CALSCR
C     ====<above step (2) and (3) should be modified>====
C     --- calculate total energy ---
C     ---- energy associated with I-th atom should be calculated ----
C      WRITE(*,*)  ' calling Eselect for dEtotal/dzi'
      CALL ESELCT(E1)
      GZ=(E1-E0)/DSTEP
C
C --- recovery of position ---
      X(I)=XMEMO
      Y(I)=YMEMO
      Z(I)=ZMEMO
C     --- Rij and Screening Factor are changed ---
C     --- Before you calculate energy next, ---
C     --- CALRIJ(step2) and CALSCR(step3) have to be called.  ---
C
C
      RETURN
      END
C######################################################################
      SUBROUTINE GO2EOF(N)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C goto end of file
C 
C There is GREAT DIFFERENT between Unix FORTRAN and ABSoft FORTRAN.
C It would be due to end line code.
C You have to write a compatible code.
C
C 1.
C END=... and ERR=... option in WRITE,READ,OPEN,,,statements.
C Difference between Unix FORTRAN and ABSoft FORTRAN
C      ----------------------------------------------
C                    Unix FORTRAN     ABSoft FORTRAN
C      ----------------------------------------------
C       END=???          ok             neglected
C       ERR=???       neglected            ok
C      ----------------------------------------------
C
C 2.
C If you use WITE(1,*) 'This is string'
C 'This is string' will be written in file in case of ABSoft FORTRAN
C ' This is string' will be written in file in case of UNIX FORTRAN
C 
C
C
      INTEGER N
      CHARACTER*17 DUM
C
  001 READ (N,'(A)',END=003,ERR=002) DUM
      GOTO 001
  002 CONTINUE
      BACKSPACE(N)
  003 CONTINUE
C
      RETURN
      END
C######################################################################
      SUBROUTINE ARBDEV
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C  ABSTRACT
C  This routine is modified from XYZDEV.
C  Direction of deviation is dedined as follows ;
C  v1:direction from 1st atom to 2nd atom,
C  v2:direction normal to v1 and nearest to 3rd atom,
C  v3:direction normal to v1 and v2.
C  1. read parameter file (and set atomic parameters)
C  2. input material name (and find from above data sets)
C  3. read spd data file
C  4. input range of deviation normalized
C  5. calculate energy as a function of deviation
C  6. end
C
      INTEGER I,NAINPF,PSYD,DIVDEV
      CHARACTER*32 SPDFIL,FILNAM
      CHARACTER*82 ONELIN
      REAL*8 NORMX(mna),NORMY(mna),NORMZ(mna),NRMDEV,D,ETOT
      REAL*8 VX(3),VY(3),VZ(3)
C
      INTEGER TOSCRN,RHOEXP
      CHARACTER*17 MEAMPF
      COMMON /SCREEN/TOSCRN /RHOEXP/RHOEXP
      COMMON /MEAMPF/MEAMPF
C
      REAL*8 X(mna),Y(mna),Z(mna)
      CHARACTER*17 NAMEOA(39),MATNAM,MATERL(mna)
      INTEGER NOATOM,PSYDAN(mna)
      REAL*8 REQUIL(39)
      COMMON /POSITN/X,Y,Z
      COMMON /NAMEOA/NAMEOA
      COMMON /MATERL/MATERL
      COMMON /NOATOM/NOATOM
      COMMON /PSYDAN/PSYDAN
      COMMON /REQUIL/REQUIL
C
C
C --- opening message ---
      WRITE(*,*) ' Calculation for Energy v.s. Deviation'
      WRITE(*,*) ' v1:direction from 1st atom to 2nd atom'
      WRITE(*,*) ' v2:direction normal to v1 and nearest to 3rd atom'
      WRITE(*,*) ' v3:direction normal to v1 and v2.'
C --- read meamPF.* ---
      CALL READPF(NAINPF)
      WRITE(*,*) ' DATA SETS IN PARAMETER FILE ARE ...'
      WRITE(*,'(4A17)')  (NAMEOA(I),I=1,NAINPF)
C --- input material name ---
  001 WRITE(*,*) ' Input material name'
      READ (*,'(A)') MATNAM
C --- find material ---
      DO 002 I=1,NAINPF
       IF(NAMEOA(I).EQ.MATNAM) THEN
          PSYD=I
          WRITE(*,*) ' Number of selected material is ',PSYD
          GOTO 003
        ENDIF
  002 CONTINUE
        WRITE(*,*) MATNAM,' is NOT seen in ',MEAMPF
        GOTO 001
C --- input file name of scalable position data ---
  003 WRITE(*,*) ' Input File Name of scalable position data'
      READ (*,'(A)') SPDFIL
C --- open the file of scalable position data ---
      OPEN (1,FILE=SPDFIL,STATUS='OLD',ERR=003)
      READ (1,'(A)') ONELIN
C     --- check the file is scalable position data or not ---
      IF (ONELIN(1:17).NE.'Scalable position') THEN
        WRITE(*,*) ' This file is not scalable position data file.'
        GOTO 003
      ENDIF
      I=0
C     --- read scalable position data ---
  004 READ (1,'(A)',ERR=005,END=005) ONELIN
      IF (ONELIN(1:1).NE.'C') THEN
        BACKSPACE(1)
        READ (1,*,ERR=003,END=003) NORMX(I+1),NORMY(I+1),NORMZ(I+1)
        WRITE(*,'(I,1P,3E20.12)') I+1,NORMX(I+1),NORMY(I+1),NORMZ(I+1)
        I=I+1
      ELSE
        WRITE(*,*) ONELIN
      ENDIF
      GOTO 004
  005 CLOSE (1)
      NOATOM=I
      WRITE(*,*) ' Number of atoms is ',NOATOM
C --- check ---
      IF (NOATOM.LT.1) THEN
        WRITE(*,*) ' Scalable position data file is DESTROYED.'
        GOTO 999
      ENDIF
C --- set v1,v2,v3 voctors ---
C     v1 can be set easily.
      VX(1)=NORMX(2)-NORMX(1)
      VY(1)=NORMY(2)-NORMY(1)
      VZ(1)=NORMZ(2)-NORMZ(1)
C     temporally set v2.
      VX(2)=NORMX(3)-NORMX(1)
      VY(2)=NORMY(3)-NORMY(1)
      VZ(2)=NORMZ(3)-NORMZ(1)
C     v3 = v1 x v2
      VX(3)=VY(1)*VZ(2)-VZ(1)*VY(2)
      VY(3)=VZ(1)*VX(2)-VX(1)*VZ(2)
      VZ(3)=VX(1)*VY(2)-VY(1)*VX(2)
C     v2 = v3 x v1
      VX(2)=VY(3)*VZ(1)-VZ(3)*VY(1)
      VY(2)=VZ(3)*VX(1)-VX(3)*VZ(1)
      VZ(2)=VX(3)*VY(1)-VY(3)*VX(1)
C --- set X,Y,Z,MATERL ---
      DO 006 I=1,NOATOM
        X(I)=REQUIL(PSYD)*NORMX(I)
        Y(I)=REQUIL(PSYD)*NORMY(I)
        Z(I)=REQUIL(PSYD)*NORMZ(I)
        MATERL(I)=MATNAM
        PSYDAN(I)=PSYD
  006 CONTINUE
C --- input normalized deviation for x,y,z position ---
C     when 0.3, deviation vectors are changed from 0.7 to 1.3.
  007 WRITE(*,*) ' Input normalized deviation (ex.0.3)'
      READ (*,*) NRMDEV
      IF ((NRMDEV.LT.0.D0).OR.(NRMDEV.GT.1.D0)) THEN
        WRITE(*,*) ' It should be in the range from 0 to 1.0.'
        GOTO 007
      ENDIF
      WRITE(*,*) ' Input division of calculation (ex.20)'
      READ (*,*) DIVDEV
        IF (DIVDEV.LE.0) GOTO 999
C --- input File Name for output, open & write head ---
  008 WRITE(*,*) ' Input File Name for out put'
      READ (*,'(A)',ERR=008) FILNAM
      OPEN (1,FILE=FILNAM,STATUS='NEW',ERR=008)
      WRITE(1,'(A)') '#  MEAM Energy v.s. 123 deviatioon'
      WRITE(1,'(A)') 'Normalized displacement of an atom'
      WRITE(1,'(A)') 'Total energy  {\it E}_{tot} (eV)'
      WRITE(1,'(2A,I3)')  
     &  '# Type of screening function : '
     & ,'default =  1 (MEAM92)       :',TOSCRN
      WRITE(1,'(2A,I3)')  
     &  '# Type of rho expansion      : '
     & ,'default =  0 (MEAM92)       :',RHOEXP
      WRITE(1,'(2A,A17)') 
     &  '# Parameter File Name        : '
     & ,'default = "meamPF" (MEAM92) :',MEAMPF
      WRITE(1,'(A31,A17)') 
     &  '#Scalable Position Data File : ',SPDFIL
      WRITE(1,'(A31,A17)') 
     &  '#Name of material            : ',MATNAM
      WRITE(1,'(A31,E17.8E3)') 
     &  '#Nearest neighbor dist (m)   : ',REQUIL(PSYD)
      WRITE(1,'(A31,I17)') 
     &  '#Deviation               +/- : ',NRMDEV
      WRITE(1,'(A31,1P,3E12.3E3)') 
     &  '#  Direction  v1 (ex,ey,ez)  : ',VX(1),VY(1),VZ(1)
      WRITE(1,'(A31,1P,3E12.3E3)') 
     &  '#  Direction  v2 (ex,ey,ez)  : ',VX(2),VY(2),VZ(2)
      WRITE(1,'(A31,1P,3E12.3E3)') 
     &  '#  Direction  v3 (ex,ey,ez)  : ',VX(3),VY(3),VZ(3)
C
      WRITE(1,'(A)') 
C
C --- roop for X start ---
      WRITE(*,'(A)') ' v1 - roop '
      WRITE(1,'(A)')       
     &'#  v1 - deviation'
      WRITE(1,'(A)')       
     &'#  delta v1/ R_i^0, total energy (eV), delta v1(m)'
      DO 100 I=-DIVDEV,DIVDEV
        D=1.0D0*NRMDEV*I/DIVDEV
        X(1)=REQUIL(PSYD)*(NORMX(1)+D*VX(1))
        Y(1)=REQUIL(PSYD)*(NORMY(1)+D*VY(1))
        Z(1)=REQUIL(PSYD)*(NORMZ(1)+D*VZ(1))
C   --- calculate Rij ---
        CALL CALRIJ
C        WRITE(*,*) ' CALRIJ terminated '
C   --- calculate screening factor ---
        CALL CALSCR
C        WRITE(*,*) ' CALSCR terminated '
C   --- calculate Etotal ---
        CALL ETOTAL(ETOT)
C        WRITE(*,*) ' ETOTAL terminated '
C   --- output to file ---
        WRITE(1,'(1P,3E18.9E3)') D,ETOT,REQUIL(PSYD)*D
  100 CONTINUE
C --- roop for X end ---
C
      WRITE(1,'(A)') 
C
C --- roop for Y start ---
      WRITE(*,'(A)') ' v2 - roop '
      WRITE(1,'(A)')       
     &'#  v2- deviation'
      WRITE(1,'(A)')       
     &'#  delta v2/ R_i^0, total energy (eV), delta v2(m)'
      DO 200 I=-DIVDEV,DIVDEV
        D=1.0D0*NRMDEV*I/DIVDEV
        X(1)=REQUIL(PSYD)*(NORMX(1)+D*VX(2))
        Y(1)=REQUIL(PSYD)*(NORMY(1)+D*VY(2))
        Z(1)=REQUIL(PSYD)*(NORMZ(1)+D*VZ(2))
C   --- calculate Rij ---
        CALL CALRIJ
C        WRITE(*,*) ' CALRIJ terminated '
C   --- calculate screening factor ---
        CALL CALSCR
C        WRITE(*,*) ' CALSCR terminated '
C   --- calculate Etotal ---
        CALL ETOTAL(ETOT)
C        WRITE(*,*) ' ETOTAL terminated '
C   --- output to file ---
        WRITE(1,'(1P,3E18.9E3)') D,ETOT,REQUIL(PSYD)*D
  200 CONTINUE
C --- roop for Y end ---
C
      WRITE(1,'(A)') 
C
C --- roop for Z start ---
      WRITE(*,'(A)') ' v3 - roop '
      WRITE(1,'(A)')       
     &'#  v3 - deviation'
      WRITE(1,'(A)')       
     &'#  delta v3/ R_i^0, total energy (eV), delta v3(m)'
      DO 300 I=-DIVDEV,DIVDEV
        D=1.0D0*NRMDEV*I/DIVDEV
        X(1)=REQUIL(PSYD)*(NORMX(1)+D*VX(3))
        Y(1)=REQUIL(PSYD)*(NORMY(1)+D*VY(3))
        Z(1)=REQUIL(PSYD)*(NORMZ(1)+D*VZ(3))
C   --- calculate Rij ---
        CALL CALRIJ
C        WRITE(*,*) ' CALRIJ terminated '
C   --- calculate screening factor ---
        CALL CALSCR
C        WRITE(*,*) ' CALSCR terminated '
C   --- calculate Etotal ---
        CALL ETOTAL(ETOT)
C        WRITE(*,*) ' ETOTAL terminated '
C   --- output to file ---
        WRITE(1,'(1P,3E18.9E3)') D,ETOT,REQUIL(PSYD)*D
  300 CONTINUE
C --- roop for Z end ---
C
      CLOSE(1)
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE XYZDEV
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C  ABSTRACT
C  1. read parameter file (and set atomic parameters)
C  2. input material name (and find from above data sets)
C  3. read spd data file
C  4. input range of deviation normalized
C  5. calculate energy as a function of deviation
C  6. end
C
      INTEGER I,NAINPF,PSYD,DIVDEV
      CHARACTER*32 SPDFIL,FILNAM
      CHARACTER*82 ONELIN
      REAL*8 NORMX(mna),NORMY(mna),NORMZ(mna),NRMDEV,D,ETOT
C
      INTEGER TOSCRN,RHOEXP
      CHARACTER*17 MEAMPF
      COMMON /SCREEN/TOSCRN /RHOEXP/RHOEXP
      COMMON /MEAMPF/MEAMPF
C
      REAL*8 X(mna),Y(mna),Z(mna)
      CHARACTER*17 NAMEOA(39),MATNAM,MATERL(mna)
      INTEGER NOATOM,PSYDAN(mna)
      REAL*8 REQUIL(39)
      COMMON /POSITN/X,Y,Z
      COMMON /NAMEOA/NAMEOA
      COMMON /MATERL/MATERL
      COMMON /NOATOM/NOATOM
      COMMON /PSYDAN/PSYDAN
      COMMON /REQUIL/REQUIL
C
C
C --- read meamPF.* ---
      CALL READPF(NAINPF)
      WRITE(*,*) ' DATA SETS IN PARAMETER FILE ARE ...'
      WRITE(*,'(4A17)')  (NAMEOA(I),I=1,NAINPF)
C --- input material name ---
  001 WRITE(*,*) ' Input material name'
      READ (*,'(A)') MATNAM
C --- find material ---
      DO 002 I=1,NAINPF
       IF(NAMEOA(I).EQ.MATNAM) THEN
          PSYD=I
          WRITE(*,*) ' Number of selected material is ',PSYD
          GOTO 003
        ENDIF
  002 CONTINUE
        WRITE(*,*) MATNAM,' is NOT seen in ',MEAMPF
        GOTO 001
C --- input file name of scalable position data ---
  003 WRITE(*,*) ' Input File Name of scalable position data'
      READ (*,'(A)') SPDFIL
C --- open the file of scalable position data ---
      OPEN (1,FILE=SPDFIL,STATUS='OLD',ERR=003)
      READ (1,'(A)') ONELIN
C     --- check the file is scalable position data or not ---
      IF (ONELIN(1:17).NE.'Scalable position') THEN
        WRITE(*,*) ' This file is not scalable position data file.'
        GOTO 003
      ENDIF
      I=0
C     --- read scalable position data ---
  004 READ (1,'(A)',ERR=005,END=005) ONELIN
      IF (ONELIN(1:1).NE.'C') THEN
        BACKSPACE(1)
        READ (1,*,ERR=003,END=003) NORMX(I+1),NORMY(I+1),NORMZ(I+1)
        WRITE(*,'(I,1P,3E20.12)') I+1,NORMX(I+1),NORMY(I+1),NORMZ(I+1)
        I=I+1
      ELSE
        WRITE(*,*) ONELIN
      ENDIF
      GOTO 004
  005 CLOSE (1)
      NOATOM=I
      WRITE(*,*) ' Number of atoms is ',NOATOM
C --- check ---
      IF (NOATOM.LT.1) THEN
        WRITE(*,*) ' Scalable position data file is DESTROYED.'
        GOTO 999
      ENDIF
C --- set X,Y,Z,MATERL ---
      DO 006 I=1,NOATOM
        X(I)=REQUIL(PSYD)*NORMX(I)
        Y(I)=REQUIL(PSYD)*NORMY(I)
        Z(I)=REQUIL(PSYD)*NORMZ(I)
        MATERL(I)=MATNAM
        PSYDAN(I)=PSYD
  006 CONTINUE
C --- input normalized deviation for x,y,z position ---
C     when 0.3, x,y,z is changed from 0.7 to 1.3.
  007 WRITE(*,*) ' Input normalized deviation (ex.0.3)'
      READ (*,*) NRMDEV
      IF ((NRMDEV.LT.0.D0).OR.(NRMDEV.GT.1.D0)) THEN
        WRITE(*,*) ' It should be in the range from 0 to 1.0.'
        GOTO 007
      ENDIF
      WRITE(*,*) ' Input division of calculation (ex.20)'
      READ (*,*) DIVDEV
        IF (DIVDEV.LE.0) GOTO 999
C --- input File Name for output, open & write head ---
  008 WRITE(*,*) ' Input File Name for out put'
      READ (*,'(A)',ERR=008) FILNAM
      OPEN (1,FILE=FILNAM,STATUS='NEW',ERR=008)
      WRITE(1,'(A)') '#  MEAM Energy v.s. XYZ deviatioon'
      WRITE(1,'(A)') 'Normalized displacement of an atom'
      WRITE(1,'(A)') 'Total energy  {\it E}_{tot} (eV)'
      WRITE(1,'(2A,I3)')  
     &               '# Type of screening function : '
     & ,'default =  1 (MEAM92)       :',TOSCRN
      WRITE(1,'(2A,I3)')  
     &               '# Type of rho expansion      : '
     & ,'default =  0 (MEAM92)       :',RHOEXP
      WRITE(1,'(2A,A17)') 
     &               '# Parameter File Name        : '
     & ,'default = "meamPF" (MEAM92) :',MEAMPF
      WRITE(1,'(A31,A17)') 
     &               '#Scalable Position Data File : ',SPDFIL
      WRITE(1,'(A31,A17)') 
     &               '#Name of material            : ',MATNAM
      WRITE(1,'(A31,E17.8E3)') 
     &               '#Nearest neighbor dist (m)   : ',REQUIL(PSYD)
      WRITE(1,'(A31,I17)') 
     &               '#Deviation               +/- : ',NRMDEV
C
      WRITE(1,'(A)') 
C
C --- roop for X start ---
      WRITE(*,'(A)') ' x - roop '
      WRITE(1,'(A)')       
     &'#  x - deviation'
      WRITE(1,'(A)')       
     &'#  delta x / R_i^0, total energy (eV), delta x (m)'
      DO 100 I=-DIVDEV,DIVDEV
        D=1.0D0*NRMDEV*I/DIVDEV
        X(1)=REQUIL(PSYD)*(NORMX(1)+D)
        Y(1)=REQUIL(PSYD)*NORMY(1)
        Z(1)=REQUIL(PSYD)*NORMZ(1)
C   --- calculate Rij ---
        CALL CALRIJ
C        WRITE(*,*) ' CALRIJ terminated '
C   --- calculate screening factor ---
        CALL CALSCR
C        WRITE(*,*) ' CALSCR terminated '
C   --- calculate Etotal ---
        CALL ETOTAL(ETOT)
C        WRITE(*,*) ' ETOTAL terminated '
C   --- output to file ---
        WRITE(1,'(1P,3E18.9E3)') D,ETOT,REQUIL(PSYD)*D
  100 CONTINUE
C --- roop for X end ---
C
      WRITE(1,'(A)') 
C
C --- roop for Y start ---
      WRITE(*,'(A)') ' y - roop '
      WRITE(1,'(A)')       
     &'#  y - deviation'
      WRITE(1,'(A)')       
     &'#  delta y / R_i^0, total energy (eV), delta y (m)'
      DO 200 I=-DIVDEV,DIVDEV
        D=1.0D0*NRMDEV*I/DIVDEV
        X(1)=REQUIL(PSYD)*NORMX(1)
        Y(1)=REQUIL(PSYD)*(NORMY(1)+D)
        Z(1)=REQUIL(PSYD)*NORMZ(1)
C   --- calculate Rij ---
        CALL CALRIJ
C        WRITE(*,*) ' CALRIJ terminated '
C   --- calculate screening factor ---
        CALL CALSCR
C        WRITE(*,*) ' CALSCR terminated '
C   --- calculate Etotal ---
        CALL ETOTAL(ETOT)
C        WRITE(*,*) ' ETOTAL terminated '
C   --- output to file ---
        WRITE(1,'(1P,3E18.9E3)') D,ETOT,REQUIL(PSYD)*D
  200 CONTINUE
C --- roop for Y end ---
C
      WRITE(1,'(A)') 
C
C --- roop for Z start ---
      WRITE(*,'(A)') ' z - roop '
      WRITE(1,'(A)')       
     &'#  z - deviation'
      WRITE(1,'(A)')       
     &'#  delta z / R_i^0, total energy (eV), delta z (m)'
      DO 300 I=-DIVDEV,DIVDEV
        D=1.0D0*NRMDEV*I/DIVDEV
        X(1)=REQUIL(PSYD)*NORMX(1)
        Y(1)=REQUIL(PSYD)*(NORMY(1)+D)
        Z(1)=REQUIL(PSYD)*NORMZ(1)
C   --- calculate Rij ---
        CALL CALRIJ
C        WRITE(*,*) ' CALRIJ terminated '
C   --- calculate screening factor ---
        CALL CALSCR
C        WRITE(*,*) ' CALSCR terminated '
C   --- calculate Etotal ---
        CALL ETOTAL(ETOT)
C        WRITE(*,*) ' ETOTAL terminated '
C   --- output to file ---
        WRITE(1,'(1P,3E18.9E3)') D,ETOT,REQUIL(PSYD)*D
  300 CONTINUE
C --- roop for Z end ---
C
      CLOSE(1)
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE SLCTPF
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
      CHARACTER*17 MEAMPF
      COMMON /MEAMPF/MEAMPF
C
C --- message & input ---
  001 WRITE(*,*) ' === Change File Name of MEAM Parameter ==='
      WRITE(*,*) ' meamPF     : Baskes,Phys.Rev.B,46,2727(1992)(MEAM92)'
      WRITE(*,*) ' meamPF.BAS : converted from "meamf" by Baskes       '
      WRITE(*,*) ' meamPF.HCP : Baskes&Johnson, Modelling Simul.Mater  '
     &                                         ,'.Sci.Eng.,2(1994)147'
      WRITE(*,*) '  Input file name ....           '
      READ (*,'(A)',ERR=001) MEAMPF
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE RHOTYP
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
      INTEGER RHOEXP
      COMMON /RHOEXP/RHOEXP
C
C --- message & input ---
  001 WRITE(*,*) ' === SetUp type of Expansion of rho ==='
      WRITE(*,*) ' 0 : eq.9a shown in Baskes,Phys.Rev.B,46,2727(1992)'
     &                                                     ,'(MEAM92)'
      WRITE(*,*) ' 1 : eq.9  shown in Baskes,Phys.Rev.B,46,2727(1992)'
      WRITE(*,*) ' 2 : eq.4c shown in Baskes&Johnson, ',
     &                     'Modelling Simul.Mater.Sci.Eng.,2(1994)147'
      WRITE(*,*) '  Input number ...   (  default is 0  )            '
      READ (*,*,ERR=001) RHOEXP
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE SCRTYP
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
      INTEGER TOSCRN
      COMMON /SCREEN/TOSCRN
C
C --- message & input ---
  001 WRITE(*,*) ' === SetUp type of screening function ==='
      WRITE(*,*) '  0 : no screening                                '
      WRITE(*,*) '  1 : MEAM92 (M.I.Baskes,Phys.Rev.B,46,2727(1992))'
      WRITE(*,*) '  2 : ... not defined yet ...                     '
      WRITE(*,*) '  Input number ...   (  default is 1  )           '
      READ (*,*,ERR=001) TOSCRN
C
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE RHOVSR
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
      REAL*8 COEFST,COEFEN,COEF,RR,RHOAL
      INTEGER COEFDV,I,NAINPF,PSYD
      CHARACTER*17 MATNAM
      CHARACTER*32 FILNAM
C
      REAL*8 REQUIL(39)
      CHARACTER*17 NAMEOA(39)
      COMMON /REQUIL/REQUIL /NAMEOA/NAMEOA
C
C --- message ---
      WRITE(*,*) ' === rho plotting vs distance R ==='
C --- read meamPF ---
      CALL READPF(NAINPF)
C --- input material name ---
  004 WRITE(*,*) ' DATA SETS IN PARAMETER FILE ARE ...'
      WRITE(*,'(4A17)')  (NAMEOA(I),I=1,NAINPF)
      WRITE(*,*) ' Input material name'
      READ (*,'(A)',ERR=004) MATNAM
C --- find material ---
      DO 005 I=1,NAINPF
        IF(NAMEOA(I).EQ.MATNAM) THEN
          PSYD=I
          WRITE(*,*) ' Number of selected material is',PSYD
          GOTO 006
        ENDIF
  005 CONTINUE
      WRITE(*,*) MATNAM,' is NOT seen in ',MEAMPF
      GOTO 004
C --- input region of distance ---
  006 WRITE(*,*) ' Input Coeff._start and Coeff._end'
      READ (*,*) COEFST,COEFEN
      WRITE(*,*) ' Input Coeff_division'
      READ (*,*) COEFDV
C --- input File Name for output, open & write head ---
  007 WRITE(*,*) ' Input File Name for out put'
      READ (*,'(A)',ERR=007) FILNAM
      OPEN (1,FILE=FILNAM,STATUS='NEW',ERR=007)
      WRITE(1,'(A)') '  rho distribution used in MEAM '
      WRITE(1,'(2A)') ' Material : ',MATNAM
      WRITE(1,'(2A)') '     R(m)        R/Requil     rho^(0)',
     &                '    rho^(1)      rho^(2)      rho^(3)'
C --- roop ---
      DO 001 I=1,COEFDV
        COEF=COEFST+(COEFEN-COEFST)*I/COEFDV
        RR=COEF*REQUIL(PSYD)
        WRITE(1,'(1P,6E13.5)') RR,COEF,RHOAL(0,PSYD,RR),
     &  RHOAL(1,PSYD,RR),RHOAL(2,PSYD,RR),RHOAL(3,PSYD,RR)
  001 CONTINUE
C --- close ---
      CLOSE(1)
C
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE ARGOFF
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C
      REAL*8 RHO0C,RHOBR0
      INTEGER NAINPF,I,PSYD
C
      REAL*8 REQUIL(39),RHO0(39)
      INTEGER ZI(39)
      CHARACTER*17   NAMEOA(39)
      CHARACTER*17 MEAMPF
      COMMON /MEAMPF/MEAMPF
      COMMON /REQUIL/REQUIL /RHO0/RHO0
      COMMON /ZI/ZI
      COMMON /NAMEOA/NAMEOA

C
      WRITE(*,*) ' =================================================='
      WRITE(*,*) ' In Table III of Baskes,Phys.Rev.B,46,2727(1992),  '
      WRITE(*,*) ' no parameter for normalization of the argument of '
      WRITE(*,*) ' embedding function F(rho).        This subroutine '
      WRITE(*,*) ' calculate the equilibrium electron density (rho_0)'
      WRITE(*,*) ' in the refference lattice.  If rho_0 is different '
      WRITE(*,*) ' from the value in "meamPF", you should be edit the'
      WRITE(*,*) ' MEAM parameter file "meamPF.*".  Good Luck !      '
      WRITE(*,*) ' =================================================='
C
C
C --- read meamPF.* ---
      WRITE(*,*) ' reading "',MEAMPF,'" ........'
      CALL READPF(NAINPF)
      WRITE(*,*)  NAINPF,'data sets are found in ',MEAMPF
      WRITE(*,*) ' DATA SETS IN PARAMETER FILE ARE ...'
      WRITE(*,'(4A17)')  (NAMEOA(I),I=1,NAINPF)
C --- roop ---
      DO 001 PSYD=1,NAINPF
C       In this program, RHOBR0(REQUIL(PSYD),PSYD)/ZI(PSYD)/RHO0(PSYD)
C       is the argument (rho) of the embedding term F(rho) for 
C       equilibrium distance in refference structure.
C       This term shall be unity. i.e.
C             RHOBR0(REQUIL(PSYD),PSYD)/ZI(PSYD)/RHO0(PSYD) = 1
        RHO0C=RHOBR0(REQUIL(PSYD),PSYD)/ZI(PSYD)
        WRITE(*,'(I3,1X,A17,1P,E23.14E3,A,E23.14E3)') 
     &            PSYD,NAMEOA(PSYD),RHO0C,' =?=',RHO0(PSYD)
  001 CONTINUE
C
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE ROSEEP
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C Rose energy plotting
      REAL*8 COEFST,COEFEN,COEF,RNEAR,ENERGY,EROSE
      INTEGER NAINPF,I,PSYD,PSYD1,COEFDV
      CHARACTER*17 MATNAM
      CHARACTER*32 FILNAM
C
      REAL*8 REQUIL(39)
      CHARACTER*17 NAMEOA(39)
      CHARACTER*17 MEAMPF
      COMMON /MEAMPF/MEAMPF
      COMMON /REQUIL/REQUIL /NAMEOA/NAMEOA
C
C --- read meamPF ---
      CALL READPF(NAINPF)
C --- input material name ---
  004 WRITE(*,*) ' DATA SETS IN PARAMETER FILE ARE ...'
      WRITE(*,'(4A17)')  (NAMEOA(I),I=1,NAINPF)
      WRITE(*,*) ' Input material name'
      READ (*,'(A)',ERR=004) MATNAM
C --- find material ---
      DO 005 I=1,NAINPF
        IF(NAMEOA(I).EQ.MATNAM) THEN
          PSYD=I
          WRITE(*,*) ' Number of selected material is',PSYD
          GOTO 006
        ENDIF
  005 CONTINUE
      WRITE(*,*) MATNAM,' is NOT seen in ',MEAMPF
      GOTO 004
C --- input region & division ---
  006 WRITE(*,*) ' Input R/Requil.(start), R/Requil.(end) (ex.0.8,3)'
      READ (*,*) COEFST,COEFEN
      WRITE(*,*) ' Input Division (ex.200)'
      READ (*,*) COEFDV
C --- input File Name for output, open & write head ---
  007 WRITE(*,*) ' Input File Name for out put'
      READ (*,'(A)',ERR=007) FILNAM
      OPEN (1,FILE=FILNAM,STATUS='NEW',ERR=007)
      WRITE(1,'(A)') ' Rose energy plotting used in MEAM '
      WRITE(1,'(2A)') ' Material : ',MATNAM
      WRITE(1,'(A)') '<-Nearest neighbor dist.-> ENERGY (eV) '
      WRITE(1,'(A)') ' Real(m)      NormalizedBy   E_Rose'
      WRITE(1,'(A)') '              Equilibrium '
      WRITE(1,'(A)') '              1st Neighbor'
C --- roop ---
      DO 100 I=0,COEFDV
        COEF=COEFST+(COEFEN-COEFST)/COEFDV*I
        RNEAR=COEF*REQUIL(PSYD)
        WRITE(*,*) ' Calculating R/Requilibrium=',COEF
C   --- calculate E_Rose ---
        PSYD1=PSYD
        ENERGY=EROSE(PSYD,PSYD1,RNEAR)
        WRITE(1,'(1P,3E13.4E3)') RNEAR,COEF,ENERGY
  100 CONTINUE
C
      CLOSE(1)
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE WRIRIJ(U)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C write Rij
C U:output unit number
C
      INTEGER U,I,J
C
      REAL*8 RIJ(mna,mna)
      COMMON /RIJ/RIJ
      REAL*8 XIJALF(mna,mna,3)
      COMMON /XIJALF/XIJALF
      INTEGER NOATOM
      COMMON /NOATOM/NOATOM
C
      WRITE(*,*) ' ########### distance Rij ###########'
C
      DO 001 I=1,NOATOM
      DO 002 J=1,NOATOM
      WRITE(U,*) ' R(',I,',',J,')=',RIJ(I,J)
  002 CONTINUE
  001 CONTINUE
C
      DO 101 I=1,NOATOM
      DO 102 J=1,NOATOM
      WRITE(U,*) ' x(',I,',',J,')=',
     &            XIJALF(I,J,1),XIJALF(I,J,2),XIJALF(I,J,3)
  102 CONTINUE
  101 CONTINUE
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE WRISCR(U)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C write Sij
C U:output unit number
C
      INTEGER U,I,J
C
      REAL*8 SCRFAC(mna,mna)
      COMMON /SCRFAC/SCRFAC
      INTEGER NOATOM
      COMMON /NOATOM/NOATOM
C
      WRITE(*,*) ' ########### Screening Factot ###########'
C
      DO 001 I=1,NOATOM
      DO 002 J=1,NOATOM
      WRITE(U,*) ' S(',I,',',J,')=',SCRFAC(I,J)
  002 CONTINUE
  001 CONTINUE
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE CALCHK
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C Energetics in uniform dilation and compression
C
      REAL*8 NORMX(mna),NORMY(mna),NORMZ(mna)
      REAL*8 COEF,COEFC
      REAL*8 EI,ETOT,RNEAR,RNEARC
      REAL*8 RHOAL,RHOIL,RHOIL0,RHOBAR,RHOBR0,FI,PHI,EROSE
      REAL*8 RHOTMP
      INTEGER I,NAINPF,PSYD,J,NTH
      CHARACTER*17 MATNAM
      CHARACTER*32 SPDFIL,FILNAM
      CHARACTER*82 ONELIN
      CHARACTER*1  VORCX,VORCY,VORCZ
C
      REAL*8 X(mna),Y(mna),Z(mna),RIJ(mna,mna)
      INTEGER NOATOM,PSYDAN(mna)
      CHARACTER*17 MATERL(mna)
      COMMON /POSITN/X,Y,Z /NOATOM/NOATOM /PSYDAN/PSYDAN /MATERL/MATERL
      COMMON /RIJ/RIJ
C
      REAL*8 RHO0(39)
      COMMON /RHO0/RHO0
      INTEGER ZI(39)
      COMMON /ZI/ZI
      REAL*8 REQUIL(39)
      CHARACTER*17 NAMEOA(39)
      COMMON /NAMEOA/NAMEOA /REQUIL/REQUIL
C
      REAL*8 EMBEDE
C
C --- input file name of scalable position data ---
  001 WRITE(*,*) ' Input File Name of scalable position data'
      READ (*,'(A)') SPDFIL
C --- open the file of scalable position data ---
      OPEN (1,FILE=SPDFIL,STATUS='OLD',ERR=001)
      READ (1,'(A)') ONELIN
C     --- check the file is scalable position data or not ---
      IF (ONELIN(1:17).NE.'Scalable position') THEN
        WRITE(*,*) ' This file is not scalable position data file.'
        GOTO 001
      ENDIF
      I=0
C     --- read scalable position data ---
  002 READ (1,'(A)',ERR=003,END=003) ONELIN
      IF (ONELIN(1:1).NE.'C') THEN
        BACKSPACE(1)
        READ (1,*,ERR=003,END=003) NORMX(I+1),NORMY(I+1),NORMZ(I+1)
        WRITE(*,'(I,1P,3E20.12)') I+1,NORMX(I+1),NORMY(I+1),NORMZ(I+1)
        I=I+1
      ELSE
        WRITE(*,*) ONELIN
      ENDIF
      GOTO 002
  003 CLOSE (1)
      NOATOM=I
      WRITE(*,*) ' Number of atoms is ',NOATOM
C --- check ---
      IF (NOATOM.LT.1) THEN
        WRITE(*,*) ' Scalable position data file is DESTROYED.'
        GOTO 999
      ENDIF
C --- read meamPF.* ---
      CALL READPF(NAINPF)
      WRITE(*,*) ' DATA SETS IN PARAMETER FILE ARE ...'
      WRITE(*,'(4A17)')  (NAMEOA(I),I=1,NAINPF)
C --- input material name ---
  004 WRITE(*,*) ' Input material name'
      READ (*,'(A)') MATNAM
C --- find material ---
      DO 005 I=1,NAINPF
       IF(NAMEOA(I).EQ.MATNAM) THEN
          PSYD=I
          WRITE(*,*) ' Number of selected material is',PSYD
          GOTO 006
        ENDIF
  005 CONTINUE
        WRITE(*,*) MATNAM,' is NOT seen in ',MEAMPF
        GOTO 004
C --- select dilation direction ---
C --- X ---
  006 WRITE(*,*) ' X-direction is Variable or Constant ? (V/C)'
      READ (*,'(A)',ERR=006) VORCX
      IF     ((VORCX.EQ.'V').OR.(VORCX.EQ.'v')) THEN
        VORCX='V'
      ELSEIF ((VORCX.EQ.'C').OR.(VORCX.EQ.'c')) THEN
        VORCX='C'
      ELSE
        GOTO 006
      ENDIF
C --- Y ---
  016 WRITE(*,*) ' Y-direction is Variable or Constant ? (V/C)'
      READ (*,'(A)',ERR=016) VORCY
      IF     ((VORCY.EQ.'V').OR.(VORCY.EQ.'v')) THEN
        VORCY='V'
      ELSEIF ((VORCY.EQ.'C').OR.(VORCY.EQ.'c')) THEN
        VORCY='C'
      ELSE
        GOTO 016
      ENDIF
C --- Z ---
  026 WRITE(*,*) ' Z-direction is Variable or Constant ? (V/C)'
      READ (*,'(A)',ERR=026) VORCZ
      IF     ((VORCZ.EQ.'V').OR.(VORCZ.EQ.'v')) THEN
        VORCZ='V'
      ELSEIF ((VORCZ.EQ.'C').OR.(VORCZ.EQ.'c')) THEN
        VORCZ='C'
      ELSE
        GOTO 026
      ENDIF
C --- input coefficient variable ( region & division ) ---
      WRITE(*,*) ' >>> For direction variable <<<'
      WRITE(*,*) ' Input Coefficient'
      READ (*,*) COEF
C --- input coefficient constant ---
      COEFC=1.D0
      IF ((VORCX.EQ.'C').OR.(VORCY.EQ.'C').OR.(VORCZ.EQ.'C')) THEN
      WRITE(*,*) ' >>> For direction constant <<<'
      WRITE(*,*) ' Input Coeff_constant (ex.1)'
      READ (*,*) COEFC
      ENDIF
C --- input File Name for output, open & write head ---
  007 WRITE(*,*) ' Input File Name for out put'
      READ (*,'(A)',ERR=007) FILNAM
      OPEN (1,FILE=FILNAM,STATUS='NEW',ERR=007)
      WRITE(1,'(A)') ' Energetics by MEAM (KunioTAKAHASHI version)'
      WRITE(1,'(A31,A17)') ' Scalable Position Data File : ',SPDFIL
      WRITE(1,'(A31,A17)') ' Name of material            : ',MATNAM
      WRITE(1,'(A)')       ' Dilation direction ( Variable or Const )'
      WRITE(1,'(A31,A1)')  '            X direction      : ',VORCX
      WRITE(1,'(A31,A1)')  '            Y direction      : ',VORCY
      WRITE(1,'(A31,A1)')  '            Z direction      : ',VORCZ
C --- only one time roop ---
        WRITE(*,*) ' Calculating R/Requilibrium=',COEF
C   --- set X,Y,Z,MATERL ---
        DO 101 J=1,NOATOM
          RNEAR=COEF*REQUIL(PSYD)
          RNEARC=COEFC*REQUIL(PSYD)
          IF (VORCX.EQ.'V') X(J)=RNEAR*NORMX(J)
          IF (VORCX.EQ.'C') X(J)=RNEARC*NORMX(J)
          IF (VORCY.EQ.'V') Y(J)=RNEAR*NORMY(J)
          IF (VORCY.EQ.'C') Y(J)=RNEARC*NORMY(J)
          IF (VORCZ.EQ.'V') Z(J)=RNEAR*NORMZ(J)
          IF (VORCZ.EQ.'C') Z(J)=RNEARC*NORMZ(J)
          MATERL(J)=MATNAM
          PSYDAN(J)=PSYD
  101   CONTINUE
C   --- calculate Rij ---
        CALL CALRIJ
        CALL WRIRIJ(1)
C        WRITE(*,*) ' CALRIJ terminated '
C   --- calculate screening factor ---
        CALL CALSCR
        CALL WRISCR(1)
C        WRITE(*,*) ' CALSCR terminated '
C   --- calculate Etotal ---
        CALL ETOTAL(ETOT)
C        WRITE(*,*) ' ETOTAL terminated '
C   --- input number of atom you glace ---
  200   WRITE(*,*) ' There are ',NOATOM,' atoms.'
        WRITE(*,*) ' Input the number of atom you glace. quit=Q'
        READ (*,'(I)',ERR=900) NTH
        IF ((NTH.LE.0).OR.(NTH.GT.NOATOM)) GOTO 200
        WRITE(1,'(A)') ' ### glaceing at',NTH,'-th atom ###'
C   --- output data ---
        WRITE(1,*) ' j      rho^{a(0)} by j-th atom'
        DO 201 J=1,NOATOM
        IF (J.NE.NTH) WRITE(1,*) J,RHOAL(0,PSYDAN(J),RIJ(NTH,J))
  201   CONTINUE
        WRITE(1,*) ' j      rho^{a(1)} by j-th atom'
        DO 202 J=1,NOATOM
        IF (J.NE.NTH) WRITE(1,*) J,RHOAL(1,PSYDAN(J),RIJ(NTH,J))
  202   CONTINUE
        WRITE(1,*) ' j      rho^{a(2)} by j-th atom'
        DO 203 J=1,NOATOM
        IF (J.NE.NTH) WRITE(1,*) J,RHOAL(2,PSYDAN(J),RIJ(NTH,J))
  203   CONTINUE
        WRITE(1,*) ' j      rho^{a(3)} by j-th atom'
        DO 204 J=1,NOATOM
        IF (J.NE.NTH) WRITE(1,*) J,RHOAL(3,PSYDAN(J),RIJ(NTH,J))
  204   CONTINUE
        WRITE(1,*) ' rho^(0) =',RHOIL(0,NTH)
        WRITE(1,*) ' rho^(1) =',RHOIL(1,NTH)
        WRITE(1,*) ' rho^(2) =',RHOIL(2,NTH)
        WRITE(1,*) ' rho^(3) =',RHOIL(3,NTH)
        WRITE(1,*) ' rho~    =',RHOBAR(NTH)
        WRITE(1,*) ' rho=rho~/Zi/rho0 =',
     &              RHOBAR(NTH)/ZI(PSYDAN(NTH))/RHO0(PSYDAN(NTH))
        WRITE(1,*) ' !!! F(rho)  =',FI(NTH)
        DO 300 J=1,NOATOM
        IF (J.NE.NTH) THEN
        WRITE(1,*) ' ### i=',NTH,' j=',J,' ###'
        WRITE(1,*) ' rho^(0)0 =',RHOIL0(0,PSYDAN(NTH),RIJ(NTH,J))
        WRITE(1,*) ' rho^(1)0 =',RHOIL0(1,PSYDAN(NTH),RIJ(NTH,J))
        WRITE(1,*) ' rho^(2)0 =',RHOIL0(2,PSYDAN(NTH),RIJ(NTH,J))
        WRITE(1,*) ' rho^(3)0 =',RHOIL0(3,PSYDAN(NTH),RIJ(NTH,J))
        WRITE(1,*) ' rho~0    =',RHOBR0(RIJ(NTH,J),PSYDAN(NTH))
        RHOTMP=RHOBR0(RIJ(NTH,J),PSYDAN(NTH))
     &         /ZI(PSYDAN(NTH))/RHO0(PSYDAN(NTH))
        WRITE(1,*) ' rho=rho~0/Zi/rho0 =',RHOTMP
        WRITE(1,*) ' !!! F(rho)  =',EMBEDE(PSYDAN(NTH),RHOTMP)
        WRITE(1,*) 'E_Rose =',EROSE(PSYDAN(NTH),PSYDAN(J),RIJ(NTH,J))
        WRITE(1,*) ' i=',NTH,' j=',J,',  phi_ij=',PHI(NTH,J)
        ENDIF
  300   CONTINUE
        WRITE(1,*) ' !!! E =',EI(NTH)
C
        GOTO 200
  900   WRITE(1,*) '  ################################'
        WRITE(1,*) '  Etotal = ',ETOT,'(eV)'
        WRITE(1,*) '  ################################'
        WRITE(*,*) '  ################################'
        WRITE(*,*) '  Etotal = ',ETOT,'(eV)'
        WRITE(*,*) '  ################################'
C
      CLOSE(1)
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE ENEGET
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C Energetics in uniform dilation and compression
C
      REAL*8 NORMX(mna),NORMY(mna),NORMZ(mna)
      REAL*8 COEFST,COEFEN,COEF,COEFC
      REAL*8 E1,E2,E3,EI,ETOT,RNEAR,RNEARC
      INTEGER I,NAINPF,PSYD,COEFDV,J
      CHARACTER*17 MATNAM
      CHARACTER*32 SPDFIL,FILNAM
      CHARACTER*82 ONELIN
      CHARACTER*1  VORCX,VORCY,VORCZ
C
      REAL*8 X(mna),Y(mna),Z(mna)
      INTEGER NOATOM,PSYDAN(mna)
      CHARACTER*17 MATERL(mna)
      COMMON /POSITN/X,Y,Z /NOATOM/NOATOM /PSYDAN/PSYDAN /MATERL/MATERL
C
      REAL*8 REQUIL(39)
      CHARACTER*17 NAMEOA(39)
      COMMON /NAMEOA/NAMEOA /REQUIL/REQUIL
C
C --- input file name of scalable position data ---
  001 WRITE(*,*) ' Input File Name of scalable position data'
      READ (*,'(A)') SPDFIL
C --- open the file of scalable position data ---
      OPEN (1,FILE=SPDFIL,STATUS='OLD',ERR=001)
      READ (1,'(A)') ONELIN
C     --- check the file is scalable position data or not ---
      IF (ONELIN(1:17).NE.'Scalable position') THEN
        WRITE(*,*) ' This file is not scalable position data file.'
        GOTO 001
      ENDIF
      I=0
C     --- read scalable position data ---
  002 READ (1,'(A)',ERR=003,END=003) ONELIN
      IF (ONELIN(1:1).NE.'C') THEN
        BACKSPACE(1)
        READ (1,*,ERR=003,END=003) NORMX(I+1),NORMY(I+1),NORMZ(I+1)
        WRITE(*,'(I,1P,3E20.12)') I+1,NORMX(I+1),NORMY(I+1),NORMZ(I+1)
        I=I+1
      ELSE
        WRITE(*,*) ONELIN
      ENDIF
      GOTO 002
  003 CLOSE (1)
      NOATOM=I
      WRITE(*,*) ' Number of atoms is ',NOATOM
C --- check ---
      IF (NOATOM.LT.1) THEN
        WRITE(*,*) ' Scalable position data file is DESTROYED.'
        GOTO 999
      ENDIF
C --- read meamPF.* ---
      CALL READPF(NAINPF)
      WRITE(*,*) ' DATA SETS IN PARAMETER FILE ARE ...'
      WRITE(*,'(4A17)')  (NAMEOA(I),I=1,NAINPF)
C --- input material name ---
  004 WRITE(*,*) ' Input material name'
      READ (*,'(A)') MATNAM
C --- find material ---
      DO 005 I=1,NAINPF
       IF(NAMEOA(I).EQ.MATNAM) THEN
          PSYD=I
          WRITE(*,*) ' Number of selected material is',PSYD
          GOTO 006
        ENDIF
  005 CONTINUE
        WRITE(*,*) MATNAM,' is NOT seen in ',MEAMPF
        GOTO 004
C --- select dilation direction ---
C --- X ---
  006 WRITE(*,*) ' X-direction is Variable or Constant ? (V/C)'
      READ (*,'(A)',ERR=006) VORCX
      IF     ((VORCX.EQ.'V').OR.(VORCX.EQ.'v')) THEN
        VORCX='V'
      ELSEIF ((VORCX.EQ.'C').OR.(VORCX.EQ.'c')) THEN
        VORCX='C'
      ELSE
        GOTO 006
      ENDIF
C --- Y ---
  016 WRITE(*,*) ' Y-direction is Variable or Constant ? (V/C)'
      READ (*,'(A)',ERR=016) VORCY
      IF     ((VORCY.EQ.'V').OR.(VORCY.EQ.'v')) THEN
        VORCY='V'
      ELSEIF ((VORCY.EQ.'C').OR.(VORCY.EQ.'c')) THEN
        VORCY='C'
      ELSE
        GOTO 016
      ENDIF
C --- Z ---
  026 WRITE(*,*) ' Z-direction is Variable or Constant ? (V/C)'
      READ (*,'(A)',ERR=026) VORCZ
      IF     ((VORCZ.EQ.'V').OR.(VORCZ.EQ.'v')) THEN
        VORCZ='V'
      ELSEIF ((VORCZ.EQ.'C').OR.(VORCZ.EQ.'c')) THEN
        VORCZ='C'
      ELSE
        GOTO 026
      ENDIF
C --- input coefficient variable ( region & division ) ---
      WRITE(*,*) ' >>> For direction variable <<<'
      WRITE(*,*) ' Input Coeff_start, Coeff_end (ex.0.8,3)'
      READ (*,*) COEFST,COEFEN
      WRITE(*,*) ' Input Coeff_division (ex.200)'
      READ (*,*) COEFDV
C --- input coefficient constant ---
      COEFC=1.D0
      IF ((VORCX.EQ.'C').OR.(VORCY.EQ.'C').OR.(VORCZ.EQ.'C')) THEN
      WRITE(*,*) ' >>> For direction constant <<<'
      WRITE(*,*) ' Input Coeff_constant (ex.1)'
      READ (*,*) COEFC
      ENDIF
C --- input File Name for output, open & write head ---
  007 WRITE(*,*) ' Input File Name for out put'
      READ (*,'(A)',ERR=007) FILNAM
      OPEN (1,FILE=FILNAM,STATUS='NEW',ERR=007)
      WRITE(1,'(A)') ' Energetics by MEAM (KunioTAKAHASHI version)'
      WRITE(1,'(A31,A17)') ' Scalable Position Data File : ',SPDFIL
      WRITE(1,'(A31,A17)') ' Name of material            : ',MATNAM
      WRITE(1,'(A)')       ' Dilation direction ( Variable or Const )'
      WRITE(1,'(A31,A1)')  '            X direction      : ',VORCX
      WRITE(1,'(A31,A1)')  '            Y direction      : ',VORCY
      WRITE(1,'(A31,A1)')  '            Z direction      : ',VORCZ
      WRITE(1,'(A31,1P,E13.5)') 
     &                     ' Coeff.for constant dirction : ',COEFC
      WRITE(1,'(2A)') '<------- Distance ------->',
     &' <--------------     ENERGY (eV)     -------------->'
      WRITE(1,'(2A)') ' between        Coeff.for ',
     &' Etotal       E1(center)   E2           E3'
      WRITE(1,'(A)')  ' 1st and 2nd    variable  '
      WRITE(1,'(A)')  ' = R_12 (m)     direction '
C --- roop ---
      DO 100 I=0,COEFDV
        COEF=COEFST+(COEFEN-COEFST)/COEFDV*I
        WRITE(*,*) ' Calculating R/Requilibrium=',COEF
C   --- set X,Y,Z,MATERL ---
        DO 101 J=1,NOATOM
          RNEAR=COEF*REQUIL(PSYD)
          RNEARC=COEFC*REQUIL(PSYD)
          IF (VORCX.EQ.'V') X(J)=RNEAR*NORMX(J)
          IF (VORCX.EQ.'C') X(J)=RNEARC*NORMX(J)
          IF (VORCY.EQ.'V') Y(J)=RNEAR*NORMY(J)
          IF (VORCY.EQ.'C') Y(J)=RNEARC*NORMY(J)
          IF (VORCZ.EQ.'V') Z(J)=RNEAR*NORMZ(J)
          IF (VORCZ.EQ.'C') Z(J)=RNEARC*NORMZ(J)
          MATERL(J)=MATNAM
          PSYDAN(J)=PSYD
  101   CONTINUE
C   --- calculate Rij ---
        CALL CALRIJ
C        WRITE(*,*) ' CALRIJ terminated '
C --- calculate screening factor ---
        CALL CALSCR
C        WRITE(*,*) ' CALSCR terminated '
C   --- calculate Etotal ---
        CALL ETOTAL(ETOT)
C        WRITE(*,*) ' ETOTAL terminated '
C   --- output data ---
        IF (NOATOM.GE.3) THEN 
          E1=EI(1)
          E2=EI(2)
          E3=EI(3)
          WRITE(1,'(1P,6E13.4E3)') RNEAR,COEF,ETOT,E1,E2,E3
        ELSE IF (NOATOM.EQ.2) THEN
          E1=EI(1)
          E2=EI(2)
          WRITE(1,'(1P,6E13.4E3)') RNEAR,COEF,ETOT,E1,E2
        ELSE IF (NOATOM.EQ.1) THEN
          E1=EI(1)
          WRITE(1,'(1P,6E13.4E3)') RNEAR,COEF,ETOT,E1
        ENDIF
C   --- check ---
C        WRITE(*,*) '  ################################'
C        WRITE(*,*) '  Etotal = ',ETOT,'(eV)'
C        WRITE(*,*) '  ################################'
  100 CONTINUE
C
      CLOSE(1)
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE TEMEAM
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C Total Energy calculation by Modified Embedded Atom Method(MEAM)
C
      REAL*8 ETOT
      INTEGER FLAG
C
C --- set NOATOM,X,Y,Z,MATERL,TOSCRN---
C     !!Caution!! TOSCRN may be changed in this subroutine 
      CALL READKT
C --- calculate Rij ---
      CALL CALRIJ
C --- calculate screening factor ---
      CALL CALSCR
C --- read meamPF.* ---
      CALL SETATP(FLAG)
      IF (FLAG.EQ.1) GOTO 999
C --- calculate Etotal ---
      CALL ETOTAL(ETOT)
C
      WRITE(*,*) '  ################################'
      WRITE(*,*) '  Etotal = ',ETOT,'(eV)'
      WRITE(*,*) '  ################################'
C
C
  999 RETURN
      END
C######################################################################
      FUNCTION EMBEDE(PSYD,NRMRHO)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C embedding energy F(rho) as a function of 
C 1.psyudo atomic number           and 
C 2.normalized atomic electrom density.
C (eq.15)
C
      INTEGER PSYD
      REAL*8 EMBEDE,NRMRHO
C
      REAL*8 ESUBLM(39),AI(39)
      COMMON /ESUBLM/ESUBLM /AI/AI
C
      IF (NRMRHO.GT.1.D-250) THEN
        EMBEDE=AI(PSYD)*ESUBLM(PSYD)*NRMRHO*DLOG(NRMRHO)
      ELSE
        EMBEDE=0.D0
      ENDIF
C
  999 RETURN
      END
C######################################################################
      FUNCTION RHOAL(L,PSYD,R)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C eq.16
C L=0,1,2,3 coresponding to s,p,d,f
      REAL*8 RHOAL,R
      REAL*8 BSTAR
      INTEGER L,PSYD
C
      REAL*8 BETA0(39),BETA1(39),BETA2(39),BETA3(39),REQUIL(39)
      COMMON /BETA/BETA0,BETA1,BETA2,BETA3
      COMMON /REQUIL/REQUIL
C
      IF      (L.EQ.0) THEN
        BSTAR=BETA0(PSYD)*(R/REQUIL(PSYD)-1.D0)
      ELSE IF (L.EQ.1) THEN
        BSTAR=BETA1(PSYD)*(R/REQUIL(PSYD)-1.D0)
      ELSE IF (L.EQ.2) THEN
        BSTAR=BETA2(PSYD)*(R/REQUIL(PSYD)-1.D0)
      ELSE IF (L.EQ.3) THEN
        BSTAR=BETA3(PSYD)*(R/REQUIL(PSYD)-1.D0)
      ELSE
        WRITE(*,*) '  Warning ! orbital type is not s,p,d, nor f.'
        WRITE(*,*) '  Check your program ; subroutine RHOAL() ! '
      ENDIF
C
      RHOAL=DEXP(-BSTAR)
C
  999 RETURN
      END
C######################################################################
      FUNCTION RHOIL0(L,PSYD,R)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C eqs. 8a,8b,8c,8d when distance is R and refference structure
C s,p,d,f contribution term  of electron density at I-th atom.
C L=0,1,2,3 coresponding to s,p,d,f as well as eqs.8a,8b,8c,8d.
C Only first neigghbours are taken into account as the same manner
C with M.I.Baskes,Phys.Rev.B,46,2727(1992).
C
      REAL*8 RHOIL0,R
      REAL*8 RHOAL
C      REAL*8 XIJ(39,3),INROOT
      INTEGER L,PSYD
C      REAL*8 S,SS
C      INTEGER I,A,B,C
C
      INTEGER ZI(39)
      CHARACTER*17 REFSTR(39)
      COMMON /ZI/ZI 
      COMMON /REFSTR/REFSTR
C
C --- set xij,yij,zij according to refference structure ---
C   where xij,      yij,      zij        are corresponding to 
C         XIJ(fn,1),XIJ(fn,2),XIJ(fn,3), 
C         fn is number of first neighbour : ZI(PSYD).
C --- dia ( ZI(PSYD)=4 ) ---
C If REFSTR(PSYD) is NOT fcc,hcp,bcc,dia NOR dim, set C XIJ(,) 
C as shown bellow and calculate the rho of refference structure
C in this proguram.
C      IF      (REFSTR(PSYD).EQ.'dia') THEN
C        XIJ(1,1)=0.D0
C        XIJ(1,2)=0.D0
C        XIJ(1,3)=1.D0
C        XIJ(2,1)=DSQRT(8.D0)/3
C        XIJ(2,2)=0.D0
C        XIJ(2,3)=-1.0D/3
C        XIJ(3,1)=-DSQRT(2.D0)/3
C        XIJ(3,2)=+DSQRT(2.D0/3)
C        XIJ(3,3)=-1.0D/3
C        XIJ(4,1)=-DSQRT(2.D0)/3
C        XIJ(4,2)=-DSQRT(2.D0/3)
C        XIJ(4,3)=-1.0D/3
C      ELSE IF (REFSTR(PSYD).EQ.'fcc') THEN
C        WRITE(*,*) ' WARNING! ',REFSTR(PSYD),'is reffered.'
C        WRITE(*,*) ' This structure type is not available in '
C        WRITE(*,*) ' subroutine RHOIL0()! Edit this subroutine !'
C      ELSE
C        WRITE(*,*) ' WARNING! ',REFSTR(PSYD),'is reffered.'
C        WRITE(*,*) ' This structure type is not available in '
C        WRITE(*,*) ' subroutine RHOIL0()! Edit this subroutine !'
C      ENDIF
C
C
C --- eq.8a ---
      IF      (L.EQ.0) THEN
        RHOIL0=RHOAL(0,PSYD,R)*ZI(PSYD)
C --- eq.8b --- =0 only for symmetric struct.(ex. dia,fcc,bcc,..)
      ELSE IF (L.EQ.1) THEN
        IF ((REFSTR(PSYD).EQ.'fcc').OR.(REFSTR(PSYD).EQ.'hcp').OR.
     &      (REFSTR(PSYD).EQ.'bcc').OR.(REFSTR(PSYD).EQ.'dia')) THEN
          RHOIL0=0.D0
        ELSE IF (REFSTR(PSYD).EQ.'dim') THEN
          RHOIL0=RHOAL(L,PSYD,R)
        ELSE
          WRITE(*,*) '  WARNING! refference structure is UNKNOWN'
          WRITE(*,*) '         subroutine RHOIL0()  L=1 '
        ENDIF
C --- eq.8c ---
      ELSE IF (L.EQ.2) THEN
        IF ((REFSTR(PSYD).EQ.'fcc').OR.(REFSTR(PSYD).EQ.'hcp').OR.
     &      (REFSTR(PSYD).EQ.'bcc').OR.(REFSTR(PSYD).EQ.'dia')) THEN
          RHOIL0=0.D0
        ELSE IF (REFSTR(PSYD).EQ.'dim') THEN
          RHOIL0=RHOAL(L,PSYD,R)*DSQRT(2.D0/3)
        ELSE
          WRITE(*,*) '  WARNING! refference structure is UNKNOWN'
          WRITE(*,*) '         subroutine RHOIL0()  L=2 '
        ENDIF
C        SS=0.D0
C        DO 203 B=1,3
C        DO 202 A=1,3
C          S=0.D0
C          DO 201 I=1,ZI(PSYD)
C          S=S+XIJ(I,A)*XIJ(I,B)
C  201     CONTINUE
C        SS=SS+S*S
C  202   CONTINUE
C  203   CONTINUE
CC
C        INROOT=SS-1.D0*ZI(PSYD)**2/3
C        IF (INROOT.LT.-1.D-10) THEN
C          WRITE(*,*) ' WARNING ARGUMENT OF SQARE ROOT TOO NEGATIVE'
C          WRITE(*,*) ' subroutine RHOIL0()! L=2'
C          INROOT=0.D0
C        ELSE IF (INROOT.LT.0.D0) THEN
C          INROOT=0.D0
C        ENDIF
CC
C        RHOIL0=DSQRT(INROOT)*RHOAL(2,PSYD,R)
C --- eq.8d ---
      ELSE IF (L.EQ.3) THEN
        IF ((REFSTR(PSYD).EQ.'fcc').OR.(REFSTR(PSYD).EQ.'bcc')) THEN
          RHOIL0=0.D0
        ELSE IF (REFSTR(PSYD).EQ.'hcp') THEN
          RHOIL0=RHOAL(L,PSYD,R)/DSQRT(3.D0)
        ELSE IF (REFSTR(PSYD).EQ.'dia') THEN
          RHOIL0=RHOAL(L,PSYD,R)*DSQRT(32.D0/9)
        ELSE IF (REFSTR(PSYD).EQ.'dim') THEN
          RHOIL0=RHOAL(L,PSYD,R)
        ELSE
          WRITE(*,*) '  WARNING! refference structure is UNKNOWN'
          WRITE(*,*) '         subroutine RHOIL0()  L=3 '
        ENDIF
C
C        SS=0.D0
C        DO 304 C=1,3
C        DO 303 B=1,3
C        DO 302 A=1,3
C          S=0.D0
C          DO 301 I=1,ZI(PSYD)
C          S=S+XIJ(I,A)*XIJ(I,B)*XIJ(I,C)
C  301     CONTINUE
C        SS=SS+S*S
C  302   CONTINUE
C  303   CONTINUE
C  304   CONTINUE
CC
C        IF (SS.LT.-1.D-10) THEN
C          WRITE(*,*) ' WARNING ARGUMENT OF SQARE ROOT TOO NEGATIVE'
C          WRITE(*,*) ' subroutine RHOIL0()! L=3'
C          SS=0.D0
C        ELSE IF (SS.LT.0.D0) THEN
C          SS=0.D0
C        ENDIF
CC
C        WRITE(*,*) ' SS=',SS
C        RHOIL0=DSQRT(SS)*RHOAL(3,PSYD,R)
C
      ELSE
        WRITE(*,*) '  Warning ! orbital type is not s,p,d, nor f.'
        WRITE(*,*) '  Check your program ; subroutine RHOIL0() ! '
      ENDIF
CC --- check ---
C      WRITE(*,*) '  In subroutin RHOIL0(L,PSYD,R)'
C      WRITE(*,*) '   L=',L,' PSYD=',PSYD,' R=',R,' RHOIL0=',RHOIL0
C
  999 RETURN
      END
C######################################################################
      FUNCTION RHOBR0(R,PSYD)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C RHOBAR when distance is R and Refference structure
C for calculation of phi_ii term from Rose energy and embedding energy.
      REAL*8 RHOBR0,R,RHOIL0
      REAL*8 S
      REAL*8 RHOI0
      INTEGER PSYD
C
      REAL*8 T0(39),T1(39),T2(39),T3(39)
      INTEGER RHOEXP
      COMMON /SPDF/T0,T1,T2,T3
      COMMON /RHOEXP/RHOEXP
C
C   --- eq.9a ---
      IF (RHOEXP.EQ.0) THEN
        RHOI0=RHOIL0(0,PSYD,R)
        S=RHOI0*(1.D0+(T1(PSYD)*(RHOIL0(1,PSYD,R)/RHOI0)**2
     &                +T2(PSYD)*(RHOIL0(2,PSYD,R)/RHOI0)**2
     &                +T3(PSYD)*(RHOIL0(3,PSYD,R)/RHOI0)**2)/2)
        IF (S.LT.-1.D-10) THEN
          WRITE(*,*) ' WARNING ARGUMENT OF SQARE ROOT TOO NEGATIVE'
          WRITE(*,*) ' subroutine RHOBR0()! S=',S
          S=0.D0
        ELSE IF (S.LT.0.D0) THEN
          S=0.D0
        ENDIF
C --- eq.9 ---
      ELSE IF (RHOEXP.EQ.1) THEN
        S=T0(PSYD)*(RHOIL0(0,PSYD,R))**2
     &   +T1(PSYD)*(RHOIL0(1,PSYD,R))**2
     &   +T2(PSYD)*(RHOIL0(2,PSYD,R))**2
     &   +T3(PSYD)*(RHOIL0(3,PSYD,R))**2
C         according to Baskes' paper, t_i^(0) 1 always unique (=1)
        IF (S.LT.-1.D-10) THEN
          WRITE(*,*) ' WARNING ARGUMENT OF SQARE ROOT TOO NEGATIVE'
          WRITE(*,*) ' subroutine RHOBR0()! S=',S
          S=0.D0
        ELSE IF (S.LT.0.D0) THEN
          WRITE(*,*) ' WARNING ARGUMENT OF SQARE ROOT NEGATIVE'
          WRITE(*,*) ' subroutine RHOBR0()! S=',S
          S=0.D0
        ENDIF
        S=DSQRT(S)
C --- eq.4c used in Baskes&Johnson,"MEAM potential for HCP metal",
C     Modelling Simul.Mater.Sci.Eng.,2(1994)147-163 ---
      ELSE IF (RHOEXP.EQ.2) THEN
        RHOI0=RHOIL0(0,PSYD,R)
        S=RHOI0*DEXP((T1(PSYD)*(RHOIL0(1,PSYD,R)/RHOI0)**2
     &               +T2(PSYD)*(RHOIL0(2,PSYD,R)/RHOI0)**2
     &               +T3(PSYD)*(RHOIL0(3,PSYD,R)/RHOI0)**2)/2)
C       according to Baskes&Johnson's paper, t_i^(0) is not be seen
C --- check of RHOEXP ---
      ELSE
        WRITE(*,*) ' WARNING UNKNOWN TYPE OF RHO EXPANSION FUNC.'
        WRITE(*,*) ' subroutine RHOBAR()!'
      ENDIF
C
      RHOBR0=S
C
  999 RETURN
      END
C######################################################################
      FUNCTION RHOIL(L,I)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C eqs. 8a,8b,8c,8d
C s,p,d,f contribution term  of electron density at I-th atom.
C L=0,1,2,3 coresponding to s,p,d,f as well as eqs.8a,8b,8c,8d.
      REAL*8 RHOIL
      REAL*8 RHOAL
      REAL*8 S,SS,SSS,INROOT
      INTEGER L,I,A,B,C
C
      INTEGER NOATOM,PSYDAN(mna)
      REAL*8 RIJ(mna,mna),XIJALF(mna,mna,3),REQUIL(39)
      COMMON /NOATOM/NOATOM /PSYDAN/PSYDAN /RIJ/RIJ /XIJALF/XIJALF
      COMMON /REQUIL/REQUIL
      REAL*8 SCRFAC(mna,mna)
      COMMON /SCRFAC/SCRFAC
C
C --- critical distance in counting all atoms except for I-th atom ---
      RCRIT=REQUIL(PSYDAN(I))*10
C
C --- eq.8a ---
      IF      (L.EQ.0) THEN
        S=0.D0
        DO 001 J=1,NOATOM
C        IF ((I.NE.J).AND.(RIJ(I,J).LT.RCRIT)) THEN
        IF (((I.NE.J).AND.(RIJ(I,J).LT.RCRIT)).AND.
     &      (SCRFAC(I,J).GT.1.D-250)) THEN
          S=S+RHOAL(0,PSYDAN(J),RIJ(I,J))*SCRFAC(I,J)
        ENDIF
  001   CONTINUE
        RHOIL=S
CC       --- check ---
C        IF (S.GT.12.D0) THEN
C          WRITE(*,*) ' WARNING THOIL(0,I) =',S,' is too large'
C          WRITE(*,*) ' I=',I
C        ENDIF
C --- eq.8b ---
      ELSE IF (L.EQ.1) THEN
        SS=0.D0
        DO 102 A=1,3
        S=0.D0
        DO 101 J=1,NOATOM
C          IF ((I.NE.J).AND.(RIJ(I,J).LT.RCRIT)) THEN
          IF (((I.NE.J).AND.(RIJ(I,J).LT.RCRIT)).AND.
     &        (SCRFAC(I,J).GT.1.D-250)) THEN
          S=S+XIJALF(I,J,A)*RHOAL(1,PSYDAN(J),RIJ(I,J))*SCRFAC(I,J)
          ENDIF
  101   CONTINUE
        SS=SS+S*S
  102   CONTINUE
C
        IF (SS.LT.-1.D-10) THEN
          WRITE(*,*) ' WARNING ARGUMENT OF SQARE ROOT TOO NEGATIVE'
          WRITE(*,*) ' subroutine RHOIL()! L=1'
          SS=0.D0
        ELSE IF (SS.LT.0.D0) THEN
          SS=0.D0
        ENDIF
C
        RHOIL=DSQRT(SS)
C --- eq.8c ---
      ELSE IF (L.EQ.2) THEN
        SS=0.D0
        DO 203 B=1,3
        DO 202 A=1,3
        S=0.D0
        DO 201 J=1,NOATOM
C          IF ((I.NE.J).AND.(RIJ(I,J).LT.RCRIT)) THEN
          IF (((I.NE.J).AND.(RIJ(I,J).LT.RCRIT)).AND.
     &        (SCRFAC(I,J).GT.1.D-250)) THEN
          S=S+XIJALF(I,J,A)*XIJALF(I,J,B)*RHOAL(2,PSYDAN(J),RIJ(I,J))
     &       *SCRFAC(I,J)
          ENDIF
  201   CONTINUE
        SS=SS+S*S
  202   CONTINUE
  203   CONTINUE
        SSS=0.D0
        DO 204 J=1,NOATOM
C          IF ((I.NE.J).AND.(RIJ(I,J).LT.RCRIT)) THEN
          IF (((I.NE.J).AND.(RIJ(I,J).LT.RCRIT)).AND.
     &        (SCRFAC(I,J).GT.1.D-250)) THEN
          SSS=SSS+RHOAL(2,PSYDAN(J),RIJ(I,J))*SCRFAC(I,J)
          ENDIF
  204   CONTINUE
        INROOT=SS-SSS*SSS/3
C
        IF (INROOT.LT.-1.D-10) THEN
          WRITE(*,*) ' WARNING ARGUMENT OF SQARE ROOT TOO NEGATIVE'
          WRITE(*,*) ' subroutine RHOIL()! L=2'
          INROOT=0.D0
        ELSE IF (INROOT.LT.0.D0) THEN
          INROOT=0.D0
        ENDIF
C
        RHOIL=DSQRT(INROOT)
C --- eq.8d ---
      ELSE IF (L.EQ.3) THEN
        SS=0.D0
        DO 304 C=1,3
        DO 303 B=1,3
        DO 302 A=1,3
        S=0.D0
        DO 301 J=1,NOATOM
C          IF ((I.NE.J).AND.(RIJ(I,J).LT.RCRIT)) THEN
          IF (((I.NE.J).AND.(RIJ(I,J).LT.RCRIT)).AND.
     &        (SCRFAC(I,J).GT.1.D-250)) THEN
          S=S+XIJALF(I,J,A)*XIJALF(I,J,B)*XIJALF(I,J,C)
     &        *RHOAL(3,PSYDAN(J),RIJ(I,J))*SCRFAC(I,J)
          ENDIF
  301   CONTINUE
        SS=SS+S*S
  302   CONTINUE
  303   CONTINUE
  304   CONTINUE
C
        IF (SS.LT.-1.D-10) THEN
          WRITE(*,*) ' WARNING ARGUMENT OF SQARE ROOT TOO NEGATIVE'
          WRITE(*,*) ' subroutine RHOIL()! L=3'
          SS=0.D0
        ELSE IF (SS.LT.0.D0) THEN
          SS=0.D0
        ENDIF
C
        RHOIL=DSQRT(SS)
      ELSE
        WRITE(*,*) '  Warning ! orbital type is not s,p,d, nor f.'
        WRITE(*,*) '  Check your program ; subroutine RHOIL() ! '
      ENDIF
C
CC --- check ---
C      WRITE(*,*) '  In subroutin RHOIL(L,I)'
C      WRITE(*,*) '   L=',L,'   I=',I,' RHOIL=',RHOIL
C
  999 RETURN
      END
C######################################################################
      FUNCTION RHOBAR(I)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C eq.9a
      REAL*8 RHOBAR,RHOIL
      REAL*8 S
      REAL*8 RHOI0
      INTEGER I
C
      INTEGER PSYDAN(mna)
      INTEGER RHOEXP
      REAL*8 T0(39),T1(39),T2(39),T3(39)
      COMMON /PSYDAN/PSYDAN
      COMMON /RHOEXP/RHOEXP
      COMMON /SPDF/T0,T1,T2,T3
C
C --- eq.9a ---
      IF (RHOEXP.EQ.0) THEN
        RHOI0=RHOIL(0,I)
        S=RHOI0*(1.D0+(T1(PSYDAN(I))*(RHOIL(1,I)/RHOI0)**2
     &                +T2(PSYDAN(I))*(RHOIL(2,I)/RHOI0)**2
     &                +T3(PSYDAN(I))*(RHOIL(3,I)/RHOI0)**2)/2)
C    according to Baskes' paper, t_i^(0) 1 always unique (=1)
C
        IF (S.LT.-1.D-10) THEN
          WRITE(*,*) ' WARNING ARGUMENT OF SQARE ROOT TOO NEGATIVE'
          WRITE(*,*) ' subroutine RHOBAR()!'
          S=0.D0
        ELSE IF (S.LT.0.D0) THEN
          S=0.D0
        ENDIF
C
C --- eq.9 ---
      ELSE IF (RHOEXP.EQ.1) THEN
        S=T0(PSYDAN(I))*(RHOIL(0,I))**2
     &   +T1(PSYDAN(I))*(RHOIL(1,I))**2
     &   +T2(PSYDAN(I))*(RHOIL(2,I))**2
     &   +T3(PSYDAN(I))*(RHOIL(3,I))**2
C
        IF (S.LT.-1.D-10) THEN
          WRITE(*,*) ' WARNING ARGUMENT OF SQARE ROOT TOO NEGATIVE'
          WRITE(*,*) ' subroutine RHOBAR()!'
          S=0.D0
        ELSE IF (S.LT.0.D0) THEN
          WRITE(*,*) ' WARNING ARGUMENT OF SQARE ROOT NEGATIVE'
          WRITE(*,*) ' subroutine RHOBAR()!'
          S=0.D0
        ENDIF
C
        S=DSQRT(S)
C --- eq.4c used in Baskes&Johnson,"MEAM potential for HCP metal",
C     Modelling Simul.Mater.Sci.Eng.,2(1994)147-163 ---
      ELSE IF (RHOEXP.EQ.2) THEN
        RHOI0=RHOIL(0,I)
        S=RHOI0*DEXP((T1(PSYDAN(I))*(RHOIL(1,I)/RHOI0)**2
     &               +T2(PSYDAN(I))*(RHOIL(2,I)/RHOI0)**2
     &               +T3(PSYDAN(I))*(RHOIL(3,I)/RHOI0)**2)/2)
C    according to Baskes&Johnson's paper, t_i^(0) is not be seen
C --- check of RHOEXP ---
      ELSE
        WRITE(*,*) ' WARNING UNKNOWN TYPE OF RHO EXPANSION FUNC.'
        WRITE(*,*) ' subroutine RHOBAR()!'
      ENDIF
C
C
      RHOBAR=S
C
C
CC --- check ---
C      IF (I.EQ.1) THEN
C        WRITE(1,*) '        RHOBAR = ',RHOBAR
C        WRITE(1,*) '  RHOIL 0,1,2,3= ',
C     &              RHOIL(0,I),RHOIL(1,I),RHOIL(2,I),RHOIL(3,I)
C        WRITE(1,*) ' T0123*RHOIL**2= ',
C     & T0(PSYDAN(I))*(RHOIL(0,I))**2,T1(PSYDAN(I))*(RHOIL(1,I))**2,
C     & T2(PSYDAN(I))*(RHOIL(2,I))**2,T3(PSYDAN(I))*(RHOIL(3,I))**2
C        WRITE(1,*) ' Sig T*RHOIL**2= ',S
C      ENDIF
C
C --- check ---
C      WRITE(*,*) '        RHOBAR = ',RHOBAR
C      WRITE(*,*) '  RHOIL 0,1,2,3= ',
C     &              RHOIL(0,I),RHOIL(1,I),RHOIL(2,I),RHOIL(3,I)
C      WRITE(*,*) ' T0123*RHOIL**2= ',
C     & T0(PSYDAN(I))*(RHOIL(0,I))**2,T1(PSYDAN(I))*(RHOIL(1,I))**2,
C     & T2(PSYDAN(I))*(RHOIL(2,I))**2,T3(PSYDAN(I))*(RHOIL(3,I))**2
C      WRITE(*,*) ' Sig T*RHOIL**2= ',S
C
C
  999 RETURN
      END
C######################################################################
      FUNCTION FI(I)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C calculate Fi(rho) term in eq.2
C embedding energy of I-th atom (I=0-mna)
      REAL*8 FI,EMBEDE
      REAL*8 RHO,RHOBAR
      INTEGER I
C
      REAL*8 RHO0(39)
      INTEGER ZI(39),PSYDAN(mna)
      COMMON /ZI/ZI
      COMMON /RHO0/RHO0
      COMMON /PSYDAN/PSYDAN
C
CC --- check ---
C      WRITE(*,*) '  FI() is called'
C
      RHO=RHOBAR(I)/ZI(PSYDAN(I))/RHO0(PSYDAN(I))
      FI=EMBEDE(PSYDAN(I),RHO)
C
CC --- check ---
C      IF (I.EQ.1) THEN
C        WRITE(1,*) '  RHO in FI(1) = ',RHO
C        WRITE(1,*) '         FI(1) = ',FI
C      ENDIF
C      WRITE(*,*) '  RHO in FI(I) = ',RHO
C      WRITE(*,*) '  F_',I,'= ',FI
C
  999 RETURN
      END
C######################################################################
      FUNCTION EROSE(PSYDI,PSYDJ,R)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
      REAL*8 EROSE,R,ASTAR
      INTEGER PSYDI,PSYDJ
C
      REAL*8 ESUBLM(39),ALPHA(39),REQUIL(39)
      COMMON /ESUBLM/ESUBLM /ALPHA/ALPHA
      COMMON /REQUIL/REQUIL
C
C --- Eii ---
      IF (PSYDI.EQ.PSYDJ) THEN
        ASTAR=ALPHA(PSYDI)*(R/REQUIL(PSYDI)-1.D0)
        EROSE=-ESUBLM(PSYDI)*(1.D0+ASTAR)*DEXP(-ASTAR)
CC    --- check ---
C        WRITE(*,*) '  EROSE called '
C        WRITE(*,*) '  R,Re          =',R,REQUIL(PSYDI)
C        WRITE(*,*) '  Alpha         =',ALPHA(PSYDI)
C        WRITE(*,*) '  a*            =',ASTAR
C        WRITE(*,*) '  exp(-a*)      =',DEXP(-ASTAR)
C        WRITE(*,*) '  (1+a*)        =',1.D0+ASTAR
C        WRITE(*,*) '  (1+a*)exp(-a*)=',(1.D0+ASTAR)*DEXP(-ASTAR)
C        WRITE(*,*) '  Eo            =',ESUBLM(PSYDI)
C        WRITE(*,*) '  E_Rose        =',EROSE
C --- Eij ---
      ELSE
        EROSE=0.D0
C
      ENDIF
CC --- check ---
C      WRITE(*,*) ' E_Rose = ',EROSE
C
  999 RETURN
      END
C######################################################################
      FUNCTION PHIII(I,J)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C eq.4 
C only first neighbours are taken into account.
      REAL*8 PHIII,EROSE,RHOBR0
      REAL*8 RHO,EMBEDE
      INTEGER I,J
C
      REAL*8 RHO0(39),RIJ(mna,mna)
      INTEGER ZI(39),PSYDAN(mna)
      COMMON /RHO0/RHO0 /RIJ/RIJ
      COMMON /ZI/ZI /PSYDAN/PSYDAN
C
C --- using rho of refference structure... ---
      RHO=RHOBR0(RIJ(I,J),PSYDAN(I))/ZI(PSYDAN(I))/RHO0(PSYDAN(I))
      PHIII=(EROSE(PSYDAN(I),PSYDAN(J),RIJ(I,J))-EMBEDE(PSYDAN(I),RHO))
     &     *2/ZI(PSYDAN(I))
CC --- check ---
C      WRITE(*,*) ' I,J in PHIII(I,J) = ',I,J
C      WRITE(*,*) ' rho in PHIII(I,J) = ',RHO,' PHIII=',PHIII
C
  999 RETURN
      END
C######################################################################
      FUNCTION PHIIJ(I,J)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
      REAL*8 PHIIJ
      INTEGER I,J
C
      PHIIJ=0.D0
C
  999 RETURN
      END
C######################################################################
      FUNCTION PHI(I,J)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
      REAL*8 PHI
      INTEGER I,J
C
      REAL*8 SCRFAC(mna,mna)
      COMMON /SCRFAC/SCRFAC
      INTEGER PSYDAN(mna)
      COMMON /PSYDAN/PSYDAN
C
      REAL*8 PHIII,PHIIJ
C
C --- case1 : screening function = 0 ---
      IF (SCRFAC(I,J).LT.1.D-250) THEN
        PHI=0.D0
C --- case2 : screening function > 0 ---
      ELSE
C      --- mono-atomic material ---
        IF (PSYDAN(I).EQ.PSYDAN(J)) THEN
          PHI=PHIII(I,J)*SCRFAC(I,J)
C      --- multi-atomic material ---
        ELSE
          PHI=PHIIJ(I,J)*SCRFAC(I,J)
        ENDIF
      ENDIF
C --- case end ---
C
  999 RETURN
      END
C######################################################################
      FUNCTION EI(I)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C
      REAL*8 EI,FI,PHI,S1,S2
      INTEGER I,J
C
      INTEGER NOATOM
      COMMON /NOATOM/NOATOM
C
CC --- check ---
C      WRITE(*,*) ' function EI() is called'
C --- embedding term ---
      S1=FI(I)
C      WRITE(*,*) ' FI() is calculated in EI()'
C --- pair-potential term ---
      S2=0.D0
      DO 002 J=1,NOATOM
      IF (I.NE.J) S2=S2+PHI(I,J)
  002 CONTINUE
      S2=S2/2
C      WRITE(*,*) ' phi-term is calculated in EI()'
C --- summation ---
      EI=S1+S2
CC --- check ---
C      WRITE(*,*) '  E_',I,'= ',EI
C      WRITE(*,*) '  F_',I,'= ',S1,'  SIGMA phi_',I,'= ',S2
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE ETOTAL(ETOT)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C  calculate total energy
C  please see "How to use..."
C
      REAL*8 ETOT,EI
      INTEGER I
C
      INTEGER NOATOM
      COMMON /NOATOM/NOATOM
C
      ETOT=0.D0
      DO 001 I=1,NOATOM
      ETOT=ETOT+EI(I)
  001 CONTINUE
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE ESELCT(ESEL)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C  calculate total energy of "selected atoms"
C  usage;
C  please replace ETOTAL to ESELCT.
C  Before calling this subroutine,
C  variables have to be set, as the case of ETOTAL
C  (step (1), (2), (3), and (4)).
C
C  When SEL(i)=0, energy of i-th atom ( EI(i) ) is NOT added to ESEL
C  When SEL(i)=1, energy of i-th atom ( EI(i) ) is     added to ESEL.
C
      REAL*8 ESEL,EI
      INTEGER I
C
      INTEGER SEL(mna)
      COMMON /SEL/SEL
      INTEGER NOATOM
      COMMON /NOATOM/NOATOM
C
      ESEL=0.D0
      DO 001 I=1,NOATOM
        IF(SEL(I).NE.0) ESEL=ESEL+EI(I)
  001 CONTINUE
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE READPF(NAINPF)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C This subroutine 
C reads MEAM parameter file defined by MEAMPF (ex.="meam.PF"),
C returns NAINPF (number of atom input from the parameter file).
C
      INTEGER I,NAINPF
      CHARACTER ONELIN*32
C
      CHARACTER*17 NAMEOA(39),REFSTR(39)
      INTEGER ZI(39),ATMNUM(39)
      REAL*8 ATMWEI(39),ALPHA(39)
      REAL*8 BETA0(39),BETA1(39),BETA2(39),BETA3(39)
      REAL*8 ESUBLM(39),AI(39),T0(39),T1(39),T2(39),T3(39)
      REAL*8 RHO0(39),REQUIL(39)
      COMMON /NAMEOA/NAMEOA
      COMMON /REFSTR/REFSTR
      COMMON /ZI/ZI
      COMMON /ATMNUM/ATMNUM
      COMMON /ATMWEI/ATMWEI
      COMMON /ALPHA/ALPHA
      COMMON /BETA/BETA0,BETA1,BETA2,BETA3
      COMMON /ESUBLM/ESUBLM
      COMMON /AI/AI
      COMMON /SPDF/T0,T1,T2,T3
      COMMON /RHO0/RHO0
      COMMON /REQUIL/REQUIL
C
      INTEGER NOATOM
      COMMON /NOATOM/NOATOM
      INTEGER PSYDAN(mna)
      COMMON /PSYDAN/PSYDAN
      CHARACTER*17 MATERL(mna)
      COMMON /MATERL/MATERL
      CHARACTER*17 MEAMPF
      COMMON /MEAMPF/MEAMPF
C
C --- read meamPF.* ---
      I=0
      OPEN(1,FILE=MEAMPF,STATUS='OLD')
  001 READ(1,'(A)',END=002,ERR=002) ONELIN
      IF ((ONELIN(1:1).EQ.'#').OR.(ONELIN(1:1).EQ.'[').OR.
     &    (ONELIN(1:1).EQ.'%').OR.(ONELIN.EQ.'')) GOTO 001
      BACKSPACE(1)
C     --- check of maxnumber start ---
      IF (I.GT.39) THEN
        WRITE(*,*) ' ######### WARNING ########'
        WRITE(*,*) ' Maximum Number of the kind of Atoms is now 39.'
        WRITE(*,*) ' However, the parameter file ',MEAMPF
        WRITE(*,*) ' has data more than 39 atoms.'
        WRITE(*,*) ' You have to modfy program.'
        STOP
      ENDIF
C     --- check of maxnumber  end  ---
      I=I+1
      READ(1,'(2A17,2I17)') NAMEOA(I),REFSTR(I),ZI(I),ATMNUM(I)
      READ(1,'(4E17.0)') ATMWEI(I),ALPHA(I),BETA0(I),BETA1(I)
      READ(1,'(4E17.0)') BETA2(I),BETA3(I),ESUBLM(I),AI(I)
      READ(1,'(4E17.0)') T0(I),T1(I),T2(I),T3(I)
      READ(1,'(2E17.0)') RHO0(I),REQUIL(I)
C
C --- check for meamPF.* reading ---
C      WRITE(*,*) ' meam parameter read from file :',MEAMPF
C      WRITE(*,'(2A17,2I17)') NAMEOA(I),REFSTR(I),ZI(I),ATMNUM(I)
C      WRITE(*,'(1P,4E17.5)') ATMWEI(I),ALPHA(I),BETA0(I),BETA1(I)
C      WRITE(*,'(1P,4E17.5)') BETA2(I),BETA3(I),ESUBLM(I),AI(I)
C      WRITE(*,'(1P,4E17.5)') T0(I),T1(I),T2(I),T3(I)
C      WRITE(*,'(1P,2E17.5)') RHO0(I),REQUIL(I)
C
      GOTO 001
  002 NAINPF=I
      CLOSE(1)
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE SETATP(FLAG)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C  Set Atomic Parameters
C  -read parameter file "meamPF.*"
C  -list up required materials and give psydo- atomic numbers to them.
C  When successfuly done,    FLAG=0
C       there are some error FLAG=1
C
      INTEGER I,J,NAINPF,FLAG
C
      CHARACTER*17 NAMEOA(39)
      COMMON /NAMEOA/NAMEOA
C
      CHARACTER*17 MATERL(mna)
      INTEGER NOATOM,PSYDAN(mna)
      COMMON /MATERL/MATERL
      COMMON /NOATOM/NOATOM /PSYDAN/PSYDAN
C
      FLAG=0
C --- read meamPF.* ---
      CALL READPF(NAINPF)
C --- set PSYDAN() start ---
      DO 100 I=1,NOATOM
      DO 101 J=1,NAINPF
        IF (MATERL(I).EQ.NAMEOA(J)) THEN
        PSYDAN(I)=J
        GOTO 102
        ENDIF
  101 CONTINUE
C     --- check ---
      WRITE(*,*) ' There is an atom whose parameters are NOT seen in',
     &' MEAM parameter file : ',MEAMPF
      WRITE(*,*) ' ',MATERL(I)
      FLAG=1
  102 CONTINUE
C     --- check of PSYDAN() numbering ---
C      WRITE(*,*) 
C     &' The',I,'-th atom ',MATERL(I),' is numbered ',PSYDAN(I),'.'
  100 CONTINUE
C --- set PSYDAN() end ---
C
  999 RETURN
      END
C######################################################################
      FUNCTION BASK92(I,K)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C screening function used in M.I.Baskes,Phys.Rev.B,46,2727(1992) 
C ( i.e. used by MEAM92 )
C Factors of the electron density contribution by K-th atom 
C which can be seen by I-th atom
C
      REAL*8 BASK92
      INTEGER I,K
C
      REAL*8 PIP6,PIP3,A,SS
      REAL*8 S,THEJIK,SIJK
      INTEGER J
C
      REAL*8 RIJ(mna,mna),XIJALF(mna,mna,3)
      COMMON /RIJ/RIJ /XIJALF/XIJALF
      REAL*8 PI
      COMMON /PI/PI
      INTEGER NOATOM
      COMMON /NOATOM/NOATOM
C
C --- const initialize ---
      PIP6=PI/6
      PIP3=PI/3
      A=500.D0
C --- factyor initialize ---
      S=1.D0
C --- multiply screening factor of other atoms (<>I-th,K-th) ---
      DO 001 J=1,NOATOM
      IF ((J.NE.I).AND.(J.NE.K)) THEN
      IF (RIJ(I,K).GE.RIJ(I,J)) THEN
C      IF ((RIJ(J,K).GE.RIJ(I,J)).AND.(RIJ(I,K).GE.RIJ(I,J))) THEN
C       --- inner product (Rij/|Rij|,Rik/|Rik|) = cos(THEJIK) ---
        THEJIK=XIJALF(I,J,1)*XIJALF(I,K,1)
     &        +XIJALF(I,J,2)*XIJALF(I,K,2)
     &        +XIJALF(I,J,3)*XIJALF(I,K,3)
C       --- ACOS(THEJIK) ---
C        WRITE(*,*) ' InnerProduct =',THEJIK
        IF     (THEJIK.GE.1.D0) THEN
          THEJIK=0.D0
        ELSEIF (THEJIK.LE.-1.D0) THEN
          THEJIK=PI
        ELSE
          THEJIK=DACOS(THEJIK)
        ENDIF
C        WRITE(*,*) ' THETA(',J,I,K,')=DACOS(InnerProduct)=',THEJIK
C       --------------------
        IF     (THEJIK.LE.PIP6) THEN
          SIJK=0.D0
        ELSEIF (THEJIK.GE.PIP3) THEN
          SIJK=1.D0
        ELSE
          SS=DSIN(PIP3-THEJIK)/DSIN(THEJIK-PIP6)
          SIJK=DEXP(-A*SS*SS*(RIJ(I,K)/RIJ(I,J)-1.D0))
        ENDIF
        S=S*SIJK
      ENDIF
      ENDIF
  001 CONTINUE
C --- fine ---
      IF (S.LT.1.D-250) S=0.D0
      BASK92=S
C      WRITE(*,*) ' BASK92(',I,K,')=',S
C
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE CALSCR
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C calculation of screening factor 
C this subroutine shall be called just after CALRIJ
C        INTEGER TOSCRN            : Type of screening function 
C              =0 : no screening
C              =1 : MEAM92(M.I.Baskes,Phys.Rev.B,46,2727(1992))
C              =2 : ... not defined yet ...
C COMMON /SCRFAC/SCRFAC
C         REAL*8  SCRFAC(mna,mna)     : Screening factor
C               defined in subroutine CALSCR
C               I-th atom sees the contribution of the electron density 
C               by J-th atom. However, The contribution is shaded by 
C               other atoms. Then, the electron density which I-th atom sees
C               is SCRFAC(I,J) x {the contiribution by J-th atom}.
C
      REAL*8 BASK92
      INTEGER I,J
C
      REAL*8 SCRFAC(mna,mna)
      INTEGER TOSCRN,NOATOM
      COMMON /SCREEN/TOSCRN /SCRFAC/SCRFAC /NOATOM/NOATOM
C
      DO 001 I=1,NOATOM
      DO 002 J=1,NOATOM
C       --- no screening ---
        IF     (TOSCRN.EQ.0) THEN
          SCRFAC(I,J)=1.D0
C       --- MEAM92(M.I.Baskes,Phys.Rev.B,46,2727(1992)) ---
        ELSEIF (TOSCRN.EQ.1) THEN
          SCRFAC(I,J)=BASK92(I,J)
C          WRITE(*,*) ' Screening at',I,' by',J,' is ',SCRFAC(I,J)
C       --- check of undefined screening ---
        ELSE
          WRITE(*,*) ' WARNING!! Number of scrennning function type'
          WRITE(*,*) '           =',TOSCRN,' is not defined !'
          WRITE(*,*) '  in subroutine CALSCR'
        ENDIF
  002 CONTINUE
  001 CONTINUE
C
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE CALRIJ
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C calculation of Rij and Rij^alpha
C XIJALF(I,J,A) = Xij^alpha. alpha=1,2,3 correspond to x,y,z.
C
C by this subroutine, following variables are set.
C COMMON /RIJ/RIJ(i,j)        :REAL*8:distance 
C                              between i-th atom and j-th atom
C COMMON /XIJALF/XIJALF(i,j,k):REAL*8:cosine to x(k=1),y(k=2),z(k=3) 
C                              of vector from i-th atom to j-th atom
C
C
C
      REAL*8 RIJVX
      REAL*8 RIJVY
      REAL*8 RIJVZ
      INTEGER I,J
C
      REAL*8 XIJALF(mna,mna,3)
      COMMON /XIJALF/XIJALF
      REAL*8 X(mna),Y(mna),Z(mna)
      COMMON /POSITN/X,Y,Z
      REAL*8 RIJ(mna,mna)
      COMMON /RIJ/RIJ
      INTEGER NOATOM
      COMMON /NOATOM/NOATOM
C
CC--- old routine (see Note "EAM and MEAM" p.21) ---
C      DO 001 I=1,NOATOM
C      DO 002 J=1,NOATOM
C        IF (I.NE.J) THEN
C        RIJVX=X(J)-X(I)
C        RIJVY=Y(J)-Y(I)
C        RIJVZ=Z(J)-Z(I)
C        RIJ(I,J)=DSQRT(RIJVX**2+RIJVY**2+RIJVZ**2)
C        XIJALF(I,J,1)=RIJVX/RIJ(I,J)
C        XIJALF(I,J,2)=RIJVY/RIJ(I,J)
C        XIJALF(I,J,3)=RIJVZ/RIJ(I,J)
C        ENDIF
C  002 CONTINUE
C  001 CONTINUE
CC--- new routine (see Note "EAM and MEAM" p.21) ---
C      WRITE(*,*) ' triangle region'
      DO 001 I=1,NOATOM-1
      DO 002 J=I+1,NOATOM
C        WRITE(*,*) ' cal RIJ vector'
        RIJVX=X(J)-X(I)
        RIJVY=Y(J)-Y(I)
        RIJVZ=Z(J)-Z(I)
C       ---distance---
C        WRITE(*,*) ' cal distance'
        RIJ(I,J)=DSQRT(RIJVX**2+RIJVY**2+RIJVZ**2)
        RIJ(J,I)=RIJ(I,J)
C       ---cosine---
C        WRITE(*,*) ' cal cosine start'
        XIJALF(I,J,1)=RIJVX/RIJ(I,J)
        XIJALF(I,J,2)=RIJVY/RIJ(I,J)
        XIJALF(I,J,3)=RIJVZ/RIJ(I,J)
C        WRITE(*,*) ' cal cosine mirror'
        XIJALF(J,I,1)=-XIJALF(I,J,1)
        XIJALF(J,I,2)=-XIJALF(I,J,2)
        XIJALF(J,I,3)=-XIJALF(I,J,3)
C        WRITE(*,*) ' cal cosine end'
  002 CONTINUE
  001 CONTINUE
CC--- new routine has been checked by Kunio 1999.5.20 ---
C
C      WRITE(*,*) ' taikaku region'
      DO 100 I=1,NOATOM
        RIJ(I,I)=0.D0
        XIJALF(I,I,1)=0.D0
        XIJALF(I,I,2)=0.D0
        XIJALF(I,I,3)=0.D0
  100 CONTINUE
C
  999 RETURN
      END
C
C######################################################################
      SUBROUTINE MODRIJ(K)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C  This has been checked by Kunio 1999.6.9 ---
C
C "RE"-calculation of Rij and Rij^alpha
C
C This routine modify only Rij concerning K-th atom
C
C
C XIJALF(I,J,A) = Xij^alpha. alpha=1,2,3 correspond to x,y,z.
C
C by this subroutine, following variables are set.
C COMMON /RIJ/RIJ(i,j)        :REAL*8:distance 
C                              between i-th atom and j-th atom
C COMMON /XIJALF/XIJALF(i,j,k):REAL*8:cosine to x(k=1),y(k=2),z(k=3) 
C                              of vector from i-th atom to j-th atom
C
C
C
      REAL*8 RIJVX
      REAL*8 RIJVY
      REAL*8 RIJVZ
      INTEGER I,K
C
      REAL*8 XIJALF(mna,mna,3)
      COMMON /XIJALF/XIJALF
      REAL*8 X(mna),Y(mna),Z(mna)
      COMMON /POSITN/X,Y,Z
      REAL*8 RIJ(mna,mna)
      COMMON /RIJ/RIJ
      INTEGER NOATOM
      COMMON /NOATOM/NOATOM
C
CC--- old routine (see Note "EAM and MEAM" p.21) ---
C      DO 001 I=1,NOATOM
C      DO 002 J=1,NOATOM
C        IF (I.NE.J) THEN
C        RIJVX=X(J)-X(I)
C        RIJVY=Y(J)-Y(I)
C        RIJVZ=Z(J)-Z(I)
C        RIJ(I,J)=DSQRT(RIJVX**2+RIJVY**2+RIJVZ**2)
C        XIJALF(I,J,1)=RIJVX/RIJ(I,J)
C        XIJALF(I,J,2)=RIJVY/RIJ(I,J)
C        XIJALF(I,J,3)=RIJVZ/RIJ(I,J)
C        ENDIF
C  002 CONTINUE
C  001 CONTINUE
CCC--- new routine (see Note "EAM and MEAM" p.21) ---
CC      WRITE(*,*) ' triangle region'
C      DO 001 I=1,NOATOM-1
C      DO 002 J=I+1,NOATOM
CC        WRITE(*,*) ' cal RIJ vector'
C        RIJVX=X(J)-X(I)
C        RIJVY=Y(J)-Y(I)
C        RIJVZ=Z(J)-Z(I)
CC       ---distance---
CC        WRITE(*,*) ' cal distance'
C        RIJ(I,J)=DSQRT(RIJVX**2+RIJVY**2+RIJVZ**2)
C        RIJ(J,I)=RIJ(I,J)
CC       ---cosine---
CC        WRITE(*,*) ' cal cosine start'
C        XIJALF(I,J,1)=RIJVX/RIJ(I,J)
C        XIJALF(I,J,2)=RIJVY/RIJ(I,J)
C        XIJALF(I,J,3)=RIJVZ/RIJ(I,J)
CC        WRITE(*,*) ' cal cosine mirror'
C        XIJALF(J,I,1)=-XIJALF(I,J,1)
C        XIJALF(J,I,2)=-XIJALF(I,J,2)
C        XIJALF(J,I,3)=-XIJALF(I,J,3)
CC        WRITE(*,*) ' cal cosine end'
C  002 CONTINUE
C  001 CONTINUE
CC--- new routine has been checked by Kunio 1999.5.20 ---
C
      DO 001 I=1,NOATOM
      IF (I.NE.K) THEN
C        WRITE(*,*) ' i=',I
        RIJVX=X(K)-X(I)
        RIJVY=Y(K)-Y(I)
        RIJVZ=Z(K)-Z(I)
C        WRITE(*,*) ' cal distance'
        RIJ(I,K)=DSQRT(RIJVX**2+RIJVY**2+RIJVZ**2)
        RIJ(K,I)=RIJ(I,K)
C        WRITE(*,*) ' cal cosine start'
        XIJALF(I,K,1)=RIJVX/RIJ(I,K)
        XIJALF(I,K,2)=RIJVY/RIJ(I,K)
        XIJALF(I,K,3)=RIJVZ/RIJ(I,K)
C        WRITE(*,*) ' cal cosine mirror'
        XIJALF(K,I,1)=-XIJALF(I,K,1)
        XIJALF(K,I,2)=-XIJALF(I,K,2)
        XIJALF(K,I,3)=-XIJALF(I,K,3)
C        WRITE(*,*) ' cal cosine end'
      ENDIF
  001 CONTINUE
C
        RIJ(K,K)=0.D0
        XIJALF(K,K,1)=0.D0
        XIJALF(K,K,2)=0.D0
        XIJALF(K,K,3)=0.D0
C
CC      WRITE(*,*) ' taikaku region'
C      DO 100 I=1,NOATOM
C        RIJ(I,I)=0.D0
C        XIJALF(I,I,1)=0.D0
C        XIJALF(I,I,2)=0.D0
C        XIJALF(I,I,3)=0.D0
C  100 CONTINUE
C
  999 RETURN
      END
C
C######################################################################
      SUBROUTINE READKT
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C read KT format file
C by this subroutine, following variables are set.
C COMMON /SCREEN/TOSCRN     :INTEGER: type of screening function.
C COMMON /POSITN/X(),Y(),Z():REAL*8 : xyz position of I-th atom
C COMMON /MATERL/MATERL()   :CHAR*17: name of I-th atom
C COMMON /NOATOM/NOATOM     :INTEGER: total number of the atoms
C
      REAL*8 X(mna),Y(mna),Z(mna)
      INTEGER TOSCRN,NOATOM,I
      CHARACTER*72 FILNAM,ONELIN
      CHARACTER*17 MATERL(mna)
      COMMON /NOATOM/NOATOM
      COMMON /POSITN/X,Y,Z
      COMMON /MATERL/MATERL
      COMMON /SCREEN/TOSCRN
C
C --- read a KT fomat file ---
      NOATOM=0
  001 WRITE(*,*) '    File name to read'
      READ (*,'(A)') FILNAM
      OPEN (1,FILE=FILNAM,STATUS='OLD',ERR=001)
C     --- skip comment line ---
  002 READ (1,'(A)') ONELIN
      WRITE(*,*) ONELIN
      IF ((ONELIN(1:1).EQ.'C').OR.(ONELIN(1:1).EQ.'c').OR.
     &    (ONELIN(1:1).EQ.'#').OR.(ONELIN(1:1).EQ.'%')) GOTO 002
C     --- read TOSCRN (type of screening function) ---
      BACKSPACE(1)
C      READ (1,'(I)',ERR=009,END=009) TOSCRN
      READ (1,*,ERR=009,END=009) TOSCRN
      WRITE(*,*) ' Type of screening is ',TOSCRN
C     --- read X(),Y(),Z(),MATERL() ---
  003 READ (1,'(A)',ERR=009,END=009) ONELIN
C      WRITE(*,*) ONELIN
      IF ((ONELIN(1:1).EQ.'C').OR.(ONELIN(1:1).EQ.'c').OR.
     &    (ONELIN(1:1).EQ.'#').OR.(ONELIN(1:1).EQ.'%')) GOTO 003
      BACKSPACE(1)
C  003 READ (1,'(3E16.7E3,1X,A2)',ERR=009,END=009) 
C  003 READ (1,*,ERR=009,END=009) 
      READ (1,*,ERR=009,END=009) 
     & X(NOATOM+1),Y(NOATOM+1),Z(NOATOM+1),MATERL(NOATOM+1)
      NOATOM=NOATOM+1
C      WRITE(*,*) ' x,y,z,mat read'
      GOTO 003
C     --- close ---
  009 CLOSE(1)
C --- check for reading data file ---
      WRITE(*,*) ' This is data in ',FILNAM
      WRITE(*,*) ' Number of atom is ',NOATOM
      DO 100 I=1,NOATOM
      WRITE(*,'(1P,3E16.7E3,1X,A10)') X(I),Y(I),Z(I),MATERL(I)
  100 CONTINUE
C
  999 RETURN
      END
C
C######################################################################
      SUBROUTINE WXYZMR(N)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C Write data file for recuction calculation
C
      INTEGER N,I
C
      REAL*8 X(mna),Y(mna),Z(mna)
      INTEGER TOSCRN,NOATOM,RHOEXP
      INTEGER REFATM(mna)
      CHARACTER*17 MEAMPF
      CHARACTER*17 MATERL(mna)
      COMMON /SCREEN/TOSCRN
      COMMON /RHOEXP/RHOEXP
      COMMON /MEAMPF/MEAMPF
      COMMON /NOATOM/NOATOM
      COMMON /POSITN/X,Y,Z
      COMMON /MATERL/MATERL
      COMMON /REFATM/REFATM
      INTEGER SEL(mna)
      COMMON /SEL/SEL
C
C --- write 3 parameters ---
      WRITE(N,'(2I7,A7,A17)') TOSCRN,RHOEXP,'       ',MEAMPF
C --- write X(),Y(),Z(),MATERL() ---
      DO 001 I=1,NOATOM
      WRITE(N,'(1P,3E17.8E3,A1,A4,I8,I8)')
     &     X(I),Y(I),Z(I),' ',MATERL(I),REFATM(I),SEL(I)
  001 CONTINUE
C --- 
C
  999 RETURN
      END
C######################################################################
      SUBROUTINE RXYZMR(FILNAM)
C######################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C Read data file for recuction calculation
C by this subroutine, following variables are set.
C COMMON /SCREEN/TOSCRN     :INTEGER: type of screening function.
C COMMON /RHOEXP/RHOEXP     :INTEGER: type of rho expansion
C COMMON /MEAMPF/MEAMPF     :CHAR*17: file name of parameter file
C COMMON /POSITN/X(),Y(),Z():REAL*8 : xyz position of I-th atom
C COMMON /MATERL/MATERL()   :CHAR*17: name of I-th atom
C COMMON /REFATM/REFATM     :INTEGER: refering atom
C COMMON /NOATOM/NOATOM     :INTEGER: total number of the atoms
C COMMON /NOMATM/NOMATM     :INTEGER: number of REFATM<0
C COMMON /SEL/SEL           :INTEGER: energy counting or not in ESELCT
C
      REAL*8 X(mna),Y(mna),Z(mna)
      INTEGER TOSCRN,NOATOM,I,RHOEXP
      INTEGER REFATM(mna)
      CHARACTER*17 MEAMPF
      CHARACTER*32 FILNAM,ONELIN
      CHARACTER*17 MATERL(mna)
      COMMON /SCREEN/TOSCRN
      COMMON /RHOEXP/RHOEXP
      COMMON /MEAMPF/MEAMPF
      COMMON /NOATOM/NOATOM
      COMMON /POSITN/X,Y,Z
      COMMON /MATERL/MATERL
      COMMON /REFATM/REFATM
      INTEGER NOMATM
      COMMON /NOMATM/NOMATM
      INTEGER SEL(mna)
      COMMON /SEL/SEL
C
C --- read a XYZMR fomat file ---
C      WRITE(*,*) ' reading ',FILNAM
C     --- open file to read ---
      OPEN (1,FILE=FILNAM,STATUS='OLD')
C     --- counter set ---
  100 NOATOM=0
C     === read TOSCRN RHOEXP MEAMPF ===
C     --- skip comment line ---
  002 READ (1,'(A)') ONELIN
      IF ((ONELIN(1:1).EQ.'C').OR.(ONELIN(1:1).EQ.'c').OR.
     &    (ONELIN(1:1).EQ.'#').OR.(ONELIN(1:1).EQ.'%')) GOTO 002
C     --- read TOSCRN RHOEXP MEAMPF ---
      BACKSPACE(1)
      READ (1,*,ERR=009) TOSCRN,RHOEXP,MEAMPF
C     === read X(),Y(),Z(),MATERL() ===
C     --- skip comment line ---
  003 READ (1,'(A)',END=019,ERR=019) ONELIN
      IF     (ONELIN(1:7).EQ.'#######') THEN
        GOTO 100
      ELSEIF ((ONELIN(1:1).EQ.'C').OR.(ONELIN(1:1).EQ.'c').OR.
     &    (ONELIN(1:1).EQ.'#').OR.(ONELIN(1:1).EQ.'%')) THEN
        GOTO 003
      ENDIF
      BACKSPACE(1)
C     --- read data line ---
      READ (1,*,END=019,ERR=009) 
     &X(NOATOM+1),Y(NOATOM+1),Z(NOATOM+1),MATERL(NOATOM+1),
     &REFATM(NOATOM+1),SEL(NOATOM+1)
      NOATOM=NOATOM+1
      GOTO 003
  009 BACKSPACE(1)
      GOTO 003
  019 CONTINUE
C     --- close ---
      CLOSE(1)
C --- count NOMATM (REFATM<0) ---
      NOMATM=0
      DO 200 I=1,NOATOM
        IF (REFATM(I).LT.0) NOMATM=NOMATM+1
  200 CONTINUE
C --- check for reading data file ---
      WRITE(*,*) ' Type of screening function (s)  :',TOSCRN
      WRITE(*,*) ' Type of rho expansion      (r)  :',RHOEXP
      WRITE(*,*) ' MEAM parameter file        (p)  :',MEAMPF
      WRITE(*,*) ' Data (s,r,p/x,y,z,m,r) file     :',FILNAM
      WRITE(*,*) ' Number of atom                  :',NOATOM
      WRITE(*,*) ' Number of atom (REFATM()<0)     :',NOMATM
C      DO 300 I=1,NOATOM
C      WRITE(*,'(1P,3E16.7E3,A1,A8,I5,I2)') 
C     &X(I),Y(I),Z(I),' ',MATERL(I),REFATM(I),SEL(I)
C  300 CONTINUE
C
  999 RETURN
      END
C
C###############################################################
      SUBROUTINE SPDMAK
C###############################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C  to make .SPD files
C
      REAL*8 X,Y,Z,NDIS(9)
      INTEGER NN,COUNT,I,J,NNTH(9),NTOT(9)
      INTEGER NOFNTH,L(255),M(255),N(255),DIST,NOFEQ(255)
      CHARACTER STR*17,FILNAM*32
C
C ------ select structure ------ 
  001 WRITE(*,*) '  Select one from ;      '
      WRITE(*,*) '  SC:  simple cubic structure'
      WRITE(*,*) '  FCC: fcc structure         '
      WRITE(*,*) '  BCC: bcc structure         '
      WRITE(*,*) '  HCP: hcp structure         '
      WRITE(*,*) '  DIA: diamond structure     '
      WRITE(*,*) '  RET: returen ...  Enter command (ex.SC)'
      READ (*,'(A)',ERR=001) STR
      IF ((STR.EQ.'ret').OR.(STR.EQ.'RET')) GOTO 999
      IF ((STR.NE.'sc').AND.(STR.NE.'SC').AND.
     &    (STR.NE.'fcc').AND.(STR.NE.'FCC').AND.
     &    (STR.NE.'bcc').AND.(STR.NE.'BCC').AND.
     &    (STR.NE.'hcp').AND.(STR.NE.'HCP').AND.
     &    (STR.NE.'dia').AND.(STR.NE.'DIA')) GOTO 001
C
C
C
C --- display message like as                                ---
C ---      cent  1st  2nd  3rd  4th  5th  6th  7th  8th  9th ---
C --- dist.   0 1.00 1.50 2.50 3.50 4.50 4.60 8.50 9.50 10.5 ---
C --- num.    1   12   12   12   12   12   12   12   16   12 ---
C --- total   1   12   24   36   48   60   72   84  100  112 ---
      CALL NEIGHB(1,STR,NNTH(1),L,M,N,DIST,NOFEQ)
      NDIS(1)=DSQRT(1.D0*DIST/6)
      NTOT(1)=1+NNTH(1)
      DO 300 I=2,9
      CALL NEIGHB(I,STR,NNTH(I),L,M,N,DIST,NOFEQ)
      NDIS(I)=DSQRT(1.D0*DIST/6)
      NTOT(I)=NTOT(I-1)+NNTH(I)
  300 CONTINUE
      WRITE(*,'(A)')
     & '     cent  1st  2nd  3rd  4th  5th  6th  7th  8th  9th'
      WRITE(*,'(A,9F5.2)') 'dist.   0',(NDIS(I),I=1,9)
      WRITE(*,'(A,9I5)')   'num.    1',(NNTH(I),I=1,9)
      WRITE(*,'(A,9I5)')   'total   1',(NTOT(I),I=1,9)
C
C
C
C --- select number --- 
  002 WRITE(*,*) '  Input Total Number of atoms (>1) '
      READ (*,*,ERR=001) NN
      IF (NN.LE.1) GOTO 001
C --- input file name for output --- 
  004 WRITE(*,*) '  Input file name for output '
      READ (*,'(A)',ERR=001) FILNAM
C --- open file for output --- 
      OPEN (1,FILE=FILNAM,STATUS='NEW',ERR=004)
      WRITE(1,'(A)') 'Scalable position data file '
      WRITE(1,'(A)') 'C'
      WRITE(1,'(A,I7,A,A3,A)') 'C',NN,' atoms in ',STR,' structure'
      WRITE(1,'(2A)') 'C The nearest neighbour distance is 1',
     &                ' (normalized).'
      WRITE(1,'(A)') 'C','C','C  center atom is zero position'
      WRITE(1,'(1P,3E24.15E3)') 0.D0,0.D0,0.D0
C --- cluster size ( number of atoms ) =1 at first for center atom ---
      COUNT=1
C --- start from 1st neighbour ---
      DO 003 I=1,NN
C       --- search NN-th neighbour's Num.,Dist.,coeff.LMN &eq.Type ---
        CALL NEIGHB(I,STR,NOFNTH,L,M,N,DIST,NOFEQ)
        COUNT=COUNT+NOFNTH
C       --- check ---
        WRITE(*,'(A,I3,A,I7,A,I3,A,I7)') 
     &        '   Num. of',I,'-th neighb.:',NOFNTH,
     &        '   dist=sqrt(',DIST,'/6)    total=',COUNT
        WRITE(1,'(A,I3,A,I7,A,I3,A,I7)') 
     &        'C  Num. of',I,'-th neighb.:',NOFNTH,
     &        '   dist=sqrt(',DIST,'/6)    total=',COUNT
CC       --- check---
C        DO 900 J=1,NOFNTH
C        WRITE(*,*) L(J),M(J),N(J),'for eq.',NOFEQ(J),' dist.=',DIST
C  900  CONTINUE
C       --- calculate x,y,z ( when lack of atoms ) ---
        IF (COUNT.LT.NN) THEN
          DO 100 J=1,NOFNTH
          CALL LMNXYZ(STR,NOFEQ(J),L(J),M(J),N(J),X,Y,Z)
C          WRITE(*,'(3I10,A,I2)') L(J),M(J),N(J),' Type :',NOFEQ(J)
C          WRITE(*,'(1P,3E24.15E3)') X,Y,Z
C          WRITE(*,'(1P,E24.15E3)')  X*X+Y*Y+Z*Z
          WRITE(1,'(1P,3E24.15E3)') X,Y,Z
  100     CONTINUE
C       --- calculate x,y,z ( when enough atoms ) ---
        ELSE
          DO 200 J=1,NOFNTH-(COUNT-NN)
          CALL LMNXYZ(STR,NOFEQ(J),L(J),M(J),N(J),X,Y,Z)
C          WRITE(*,'(3I10,A,I2)') L(J),M(J),N(J),' Type :',NOFEQ(J)
C          WRITE(*,'(1P,3E24.15E3)') X,Y,Z
C          WRITE(*,'(1P,E24.15E3)')  X*X+Y*Y+Z*Z
          WRITE(1,'(1P,3E24.15E3)') X,Y,Z
  200     CONTINUE
          GOTO 990
        ENDIF
  003 CONTINUE
C
C
C
  990 CLOSE(1)
      GOTO 001
C
  999 RETURN
      END
C
C####################################################################
      SUBROUTINE LMNXYZ(STR,NNOFEQ,LL,MM,NN,XX,YY,ZZ)
C####################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C If you set STR, NNOFEQ, LL, MM, NN and call this subroutine
C (XX,YY,ZZ) = LL e1 + MM e2 + NN e3 + d
C              where e1,e2,e3 is unit vector defined for each STR
C              and d is a 0 vector NNOFEQ=1
C              and d is a NON 0 vector when NNOFEQ=2
C
      REAL*8 EX1,EY1,EZ1,EX2,EY2,EZ2,EX3,EY3,EZ3,DX,DY,DZ
      REAL*8 XX,YY,ZZ
      INTEGER LL,MM,NN
      INTEGER NNOFEQ
      CHARACTER STR*17
C
CC --- check ---
C      WRITE(*,'(3I10,A,I2,1X,A)') LL,MM,NN,' Type :',NNOFEQ,STR
C --- initialize ---
      DX=0.D0
      DY=0.D0
      DZ=0.D0
C --- set Unit Vector ---
C --- simple cubic ---
      IF      ((STR.EQ.'sc').OR.(STR.EQ.'SC')) THEN
        EX1=1.D0
        EY1=0.D0
        EZ1=0.D0
        EX2=0.D0
        EY2=1.D0
        EZ2=0.D0
        EX3=0.D0
        EY3=0.D0
        EZ3=1.D0
C       --- chek ---
          IF (NNOFEQ.NE.1) THEN
          WRITE(*,*) '  WARNING !! type of equation is only 1 for SC'
          WRITE(*,*) '    in subroutine CVE123()'
          ENDIF
C --- face center cubic ---
      ELSE IF ((STR.EQ.'fcc').OR.(STR.EQ.'FCC')) THEN
        EX1=1.D0/DSQRT(2.D0)
        EY1=1.D0/DSQRT(2.D0)
        EZ1=0.D0
        EX2=0.D0
        EY2=1.D0/DSQRT(2.D0)
        EZ2=1.D0/DSQRT(2.D0)
        EX3=1.D0/DSQRT(2.D0)
        EY3=0.D0
        EZ3=1.D0/DSQRT(2.D0)
C       --- chek ---
          IF (NNOFEQ.NE.1) THEN
          WRITE(*,*) '  WARNING !! type of equation is only 1 for FCC'
          WRITE(*,*) '    in subroutine CVE123()'
          ENDIF
C --- body center cubic ---
      ELSE IF ((STR.EQ.'bcc').OR.(STR.EQ.'BCC')) THEN
        EX1=-1.D0/DSQRT(3.D0)
        EY1=+1.D0/DSQRT(3.D0)
        EZ1=+1.D0/DSQRT(3.D0)
        EX2=+1.D0/DSQRT(3.D0)
        EY2=-1.D0/DSQRT(3.D0)
        EZ2=+1.D0/DSQRT(3.D0)
        EX3=+1.D0/DSQRT(3.D0)
        EY3=+1.D0/DSQRT(3.D0)
        EZ3=-1.D0/DSQRT(3.D0)
C       --- chek ---
          IF (NNOFEQ.NE.1) THEN
          WRITE(*,*) '  WARNING !! type of equation is only 1 for BCC'
          WRITE(*,*) '    in subroutine CVE123()'
          ENDIF
C --- hexagonal closed pack ---
      ELSE IF ((STR.EQ.'hcp').OR.(STR.EQ.'HCP')) THEN
        EX1=1.D0
        EY1=0.D0
        EZ1=0.D0
        EX2=1.D0/2
        EY2=DSQRT(3.D0)/2
        EZ2=0.D0
        EX3=0.D0
        EY3=0.D0
        EZ3=2.D0*DSQRT(2.D0/3)
C       --- chek & set d vector ---
          IF (NNOFEQ.EQ.1) THEN
            DX=0.D0
            DY=0.D0
            DZ=0.D0
          ELSE IF (NNOFEQ.EQ.2) THEN
            DX=1.D0/2.D0
            DY=1.D0/2.D0/DSQRT(3.D0)
            DZ=DSQRT(2.D0/3.D0)
          ELSE
            WRITE(*,*) '  WARNING !! UNdefined type of equation called'
            WRITE(*,*) '  for HCP  in subroutine CVE123()'
          ENDIF
C --- diamond ---
      ELSE IF ((STR.EQ.'dia').OR.(STR.EQ.'DIA')) THEN
        EX1=0.D0
        EY1=2.D0/DSQRT(3.D0)
        EZ1=2.D0/DSQRT(3.D0)
        EX2=2.D0/DSQRT(3.D0)
        EY2=0.D0
        EZ2=2.D0/DSQRT(3.D0)
        EX3=2.D0/DSQRT(3.D0)
        EY3=2.D0/DSQRT(3.D0)
        EZ3=0.D0
C       --- chek & set d vector ---
          IF (NNOFEQ.EQ.1) THEN
            DX=0.D0
            DY=0.D0
            DZ=0.D0
          ELSE IF (NNOFEQ.EQ.2) THEN
            DX=1.D0/DSQRT(3.D0)
            DY=1.D0/DSQRT(3.D0)
            DZ=1.D0/DSQRT(3.D0)
          ELSE
            WRITE(*,*) '  WARNING !! UNdefined type of equation called'
            WRITE(*,*) '  for HCP  in subroutine LMNXYZ()'
          ENDIF
C --- check of undefined structure if it were called ---
      ELSE
        WRITE(*,*) '  WARNING !! UNdefined type of lattice is called.'
        WRITE(*,*) '    in subroutine LMNXYZ()'
      ENDIF
C
C --- product l,m,n to unit vector to calculate x,y,z ---
      XX=EX1*LL+EX2*MM+EX3*NN+DX
      YY=EY1*LL+EY2*MM+EY3*NN+DY
      ZZ=EZ1*LL+EZ2*MM+EZ3*NN+DZ
C
C
      RETURN
      END
C####################################################################
      SUBROUTINE NEIGHB(NN,STR,NOFNTH,NTHL,NTHM,NTHN,DIST,NOFEQ)
C####################################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
C If you set NN & STR and call this subroutine
C NOFNTH         : Number of NN-th neighbour atom
C NTHL,NTHM,NTHN : lattice position parameter, i.e.,
C                : r = NTHL e_1 + NTHM e_2 + NTHN e_3      (NOFEQ=1)
C                : r = NTHL e_1 + NTHM e_2 + NTHN e_3 + r0 (NOFEQ=2)
C                : where e_1, e_2, e_3 is defined for each STR.
C                : There are 2 eqs. for hcp and diamond STR,
C                : which is presented with NOFEQ.
C DIST           : distance normalized by 1st neighbour distance 
C                : is sqrt(DIST/6).
C NOFEQ          : expressed above...
C are returned.
      INTEGER NN,NOFNTH,NTHL(255),NTHM(255),NTHN(255),DIST,NOFEQ(255)
      INTEGER NSEARC,L,M,N,I,NTH,FLAG
      CHARACTER STR*17
C
C --- N_search ---
      NSEARC=1
C --- counting flag : number of neighbour when >0---
      FLAG=0
C --- searching number (search NTH-th neighbour) ---
  100 NTH=0
C --- try to count by using N_search(NSEARC) ---
      DO 010 I=1,MAXCNT(NSEARC,STR)
C      WRITE(*,*) ' I=',I
C --- search L,M,N=-N_search,N_search ---
        DO 001 L=-NSEARC,NSEARC
        DO 002 M=-NSEARC,NSEARC
        DO 003 N=-NSEARC,NSEARC
          IF (INRT1(L,M,N,STR).EQ.I) THEN
C            WRITE(*,*) ' FOUND! L,M,N=',L,M,N
            IF(FLAG.EQ.0)THEN
            NTH=NTH+1
            ENDIF
            FLAG=FLAG+1
              NTHL(FLAG)=L
              NTHM(FLAG)=M
              NTHN(FLAG)=N
              DIST=I
              NOFEQ(FLAG)=1
          ENDIF
          IF ((STR.EQ.'hcp').OR.(STR.EQ.'dia').OR.
     &        (STR.EQ.'HCP').OR.(STR.EQ.'DIA')) THEN
          IF (INRT2(L,M,N,STR).EQ.I) THEN
C            WRITE(*,*) ' FOUND! L,M,N=',L,M,N
            IF(FLAG.EQ.0)THEN
            NTH=NTH+1
            ENDIF
            FLAG=FLAG+1
              NTHL(FLAG)=L
              NTHM(FLAG)=M
              NTHN(FLAG)=N
              DIST=I
              NOFEQ(FLAG)=2
          ENDIF
          ENDIF
  003   CONTINUE
  002   CONTINUE
  001   CONTINUE
C
      IF (FLAG.GE.255) THEN
        WRITE(*,*) ' WARNING ! number of N-th neighbour over flow'
        WRITE(*,*) '  variables are DESTROYED !! stop program !! '
        WRITE(*,*) '  change MAX number of NTHL(),NTHM(),NTHN()'
        WRITE(*,*) '  and re-compile !'
        WRITE(*,*) '  subroutine NEIGHB()'
      ENDIF
C
      IF (NTH.EQ.NN) THEN
        NOFNTH=FLAG
        GOTO 999
      ENDIF
C
      FLAG=0
  010 CONTINUE
      NSEARC=NSEARC+1
      GOTO 100
C
  999 RETURN
      END
C###############################################################
      FUNCTION MAXCNT(NN,STR)
C###############################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
      INTEGER MAXCNT,NN
      CHARACTER STR*17
C
      IF      ((STR.EQ.'sc').OR.(STR.EQ.'SC')) THEN
        MAXCNT=6*NN*NN-1
      ELSE IF ((STR.EQ.'fcc').OR.(STR.EQ.'FCC')) THEN
        MAXCNT=4*NN*NN-1
      ELSE IF ((STR.EQ.'hcp').OR.(STR.EQ.'HCP')) THEN
        MAXCNT=2*INT((3.0D0/2*NN-0.5D0)**2)-1
      ELSE IF ((STR.EQ.'bcc').OR.(STR.EQ.'BCC')) THEN
        MAXCNT=2*NN*NN-1
      ELSE IF ((STR.EQ.'dia').OR.(STR.EQ.'DIA')) THEN
        MAXCNT=INT(3.0D0/3*(2*NN-1)**2)-1
      ELSE
        MAXCNT=-1
      ENDIF
C
      RETURN
      END
C###############################################################
      FUNCTION INRT1(L,M,N,STR)
C###############################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
      INTEGER INRT1,L,M,N
      CHARACTER STR*17
C
      IF      ((STR.EQ.'sc').OR.(STR.EQ.'SC')) THEN
        INRT1=(L*L+M*M+N*N)*6
      ELSE IF ((STR.EQ.'fcc').OR.(STR.EQ.'FCC')) THEN
        INRT1=((L+M)**2+(M+N)**2+(N+L)**2)*3
      ELSE IF ((STR.EQ.'hcp').OR.(STR.EQ.'HCP')) THEN
        INRT1=(3*L*L+3*L*M+3*M*M+8*N*N)*2
      ELSE IF ((STR.EQ.'bcc').OR.(STR.EQ.'BCC')) THEN
        INRT1=((-L+M+N)**2+(L-M+N)**2+(L+M-N)**2)*2
      ELSE IF ((STR.EQ.'dia').OR.(STR.EQ.'DIA')) THEN
        INRT1=((L+M)**2+(M+N)**2+(N+L)**2)*8
      ELSE
        INRT1=-1
      ENDIF
C
      RETURN
      END
C###############################################################
      FUNCTION INRT2(L,M,N,STR)
C###############################################################
      parameter (mna=1023,mncmn=17,mnka=39,sr8n=1.D-250)
      INTEGER INRT2,L,M,N
      CHARACTER STR*17
C
      IF      ((STR.EQ.'hcp').OR.(STR.EQ.'HCP')) THEN
        INRT2=(3*L*L+3*L*M+3*M*M+8*N*N+3*L+3*M+8*N+3)*2
      ELSE IF ((STR.EQ.'dia').OR.(STR.EQ.'DIA')) THEN
        INRT2=((2*L+2*M+1)**2+(2*M+2*N+1)**2+(2*N+2*L+1)**2)*2
      ELSE
        INRT2=-1
      ENDIF
C
      RETURN
      END
