I think I've also fixed the gfortran issue. However, since this requires changes in x_lapw, I cannot put this into the repository.

As mentioned, there is a new version in the download area of www.wien2k.at. Under files you can download SRC_pes.tar.gz

For gfortran you should in addition replace the two attached files
periodic_table.f and estimate_cross.f

AND  change the pes: section in x_lapw to:

case pes:
set exe = pes
cat << theend > $def
2, '$WIENROOT/SRC_pes/Crosssections/', 'unknown','formatted',0
4, '$WIENROOT/SRC_pes/Crosssections/cross.txt','old','formatted',0
23, 'pes.template',      'unknown',    'formatted',0                <--
24, 'valence.reconfig',      'unknown',    'formatted',0            <--
....



On 3/13/20 2:14 PM, Peter Blaha wrote:
Daer Pavel,

Thanks for your report.

I tried to reproduce it, but my version of pes has already been changed significantly. In particular optimize_charge.f is now quite different.

However, when compiling with -C I detected another problem with the variable "start" (never set to zero) and in read_dos.f.

I've uploaded SRC_pes.tar.gz into the "files" directory of the wien2k-download and you can download it from there, if interested.

PS: I've not fixed the gfortran problem yet, but will try to do it soon. So if it is not timely, you can also wait with the download until this is also fixed.

PPS: Since you have N-s basis and it is used in the fit, you must include its PDOS also. Thus you must include the PDOS for lower energy, such that the N-s band is included. Otherwise the fit may give nonsense.

Best regards
Peter


On 3/13/20 11:58 AM, Pavel Ondračka wrote:
Dear Wien2k mailing list,

I'm experiencing a crash when trying to calculate valence band spectra
for VN.

(This is a resend of previous email which is stuck in the queue due to
being slightly over the size limit, now with a link instead. I
apologize for double posting if the original one eventually makes it to
the list as well.)

There is out of bounds write during optimization of q_spheres:
Program received signal SIGSEGV, Segmentation fault.
0x000000000040d494 in optimize_charge () at optimize_charge.f:95
95                 temp(l,recon_counter)=temp(l,j)+temp(l,i)
(gdb) print recon_counter
$1 = 27
(gdb) print output_counter
$2 = 24
(it tries to write at index 27) but the size is just 24 (defined
by output_counter)

The files needed to reproduce this and the terminal output (together
with the manual keyboard input needed to reproduce the crash) are at
https://drive.google.com/open?id=1NZ8lSkfrgigtdQZrDZLp8Y-mFf4uSyk_
. I'm not a regular user of the pes program so there is a high
chance that there is something wrong with my input.

BTW While taking a quick look I spotted some likely unrelated small
issues, for instance pes is also influenced by the well known issue
with gfortran using the units 5 and 6 (have to change it manually to
something else otherwise stdin and stdout doesn't work) and there are
some valgrind warnings even before the crash, for example:

==57304== Conditional jump or move depends on uninitialised value(s)
==57304==    at 0x419919: abs_smooth_ (SPLINE.f:173)
==57304==    by 0x41852E: setup_ (SPLINE.f:91)
==57304==    by 0x4199FE: spline_ (SPLINE.f:16)
==57304==    by 0x41582B: read_database2_ (read_database2.f:124)
==57304==    by 0x403FE7: MAIN__ (pes.f:151)
==57304==    by 0x4066FF: main (pes.f:3)
==57304==  Uninitialised value was created by a stack allocation
==57304==    at 0x4199B3: spline_ (SPLINE.f:1)

SPLINE.f:173
    if (x >= delta_x) then
The unuinitialized variable is the delta_x which was passed from setup
(SPLINE.f:91):
call abs_smooth(m4 - m3, delta_x, w1)
and was itself allocated on the stack at the beginning of spline but is
not initialized anywhere as far as I can see. So I set it to 0.0d0 (the
default for ifort).

and one more which should be harmless...
==57389== Conditional jump or move depends on uninitialised value(s)
==57389==    at 0x47213EC5: bcmp (vg_replace_strmem.c:1113)
==57389==    by 0x474A3B7A: _gfortran_compare_string
(string_intrinsics_inc.c:98)
==57389==    by 0x41677F: read_outputst_ (read_outputst.f:37)
==57389==    by 0x4048E2: MAIN__ (pes.f:222)
==57389==    by 0x4066FF: main (pes.f:3)
==57389==  Uninitialised value was created by a stack allocation
==57389==    at 0x416559: read_outputst_ (read_outputst.f:1)

Best regards
Pavel

_______________________________________________
Wien mailing list
Wien@zeus.theochem.tuwien.ac.at
http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
SEARCH the MAILING-LIST at: http://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/index.html



--

                                      P.Blaha
--------------------------------------------------------------------------
Peter BLAHA, Inst.f. Materials Chemistry, TU Vienna, A-1060 Vienna
Phone: +43-1-58801-165300             FAX: +43-1-58801-165982
Email: bl...@theochem.tuwien.ac.at    WIEN2k: http://www.wien2k.at
WWW:   http://www.imc.tuwien.ac.at/TC_Blaha
--------------------------------------------------------------------------
SUBROUTINE periodic_table

       use p_table,only : atom,composition
       use struct, only : mult, nat,aname

       implicit none
       integer     :: i,j,n,counter,lt,ios
       logical     :: exist
       character*1 :: answer,temp1,temp2
       character*2,dimension(:,:),allocatable ::reconf
       character*80 string
       atom   = '  '
!! The subroutine periodic_table reads atoms valence partial orbital names (according to provided periodic table data) into an array (composition). The periodic table data could be modified to include/exclude partial orbitals.       

       atom(1,1) ='H'  ;atom(1,2) ='1s'                                                   
       atom(2,1) ='He' ;atom(2,2) ='1s'
                                                         
       atom(3,1) ='Li' ;atom(3,2) ='2s'                                          
       atom(4,1) ='Be' ;atom(4,2) ='2s'                                                          
       atom(5,1) ='B'  ;atom(5,2) ='2s' ;atom(5,3) ='2p'                                        
       atom(6,1) ='C'  ;atom(6,2) ='2s' ;atom(6,3) ='2p'                                        
       atom(7,1) ='N'  ;atom(7,2) ='2s' ;atom(7,3) ='2p'                                         
       atom(8,1) ='O'  ;atom(8,2) ='2s' ;atom(8,3) ='2p'                                         
       atom(9,1) ='F'  ;atom(9,2) ='2s' ;atom(9,3) ='2p'                                         
       atom(10,1)='Ne' ;atom(10,2)='2s' ;atom(10,3)='2p' 
                                       
       atom(11,1)='Na' ;atom(11,2)='3s' ;atom(11,3)='3p'                       
       atom(12,1)='Mg' ;atom(12,2)='3s' ;atom(12,3)='3p'                                        
       atom(13,1)='Al' ;atom(13,2)='3s' ;atom(13,3)='3p'                                         
       atom(14,1)='Si' ;atom(14,2)='3s' ;atom(14,3)='3p'                                        
       atom(15,1)='P'  ;atom(15,2)='3s' ;atom(15,3)='3p'                                         
       atom(16,1)='S'  ;atom(16,2)='3s' ;atom(16,3)='3p'                                         
       atom(17,1)='Cl' ;atom(17,2)='3s' ;atom(17,3)='3p'                                         
       atom(18,1)='Ar' ;atom(18,2)='3s' ;atom(18,3)='3p'

       atom(19,1)='K'  ;atom(19,2)='4s' ;atom(19,3)='3p' 
       atom(20,1)='Ca' ;atom(20,2)='4s' ;atom(20,3)='4p' ;atom(20,4)='3d'       
       atom(21,1)='Sc' ;atom(21,2)='4s' ;atom(21,3)='4p' ;atom(21,4)='3d'     
       atom(22,1)='Ti' ;atom(22,2)='4s' ;atom(22,3)='4p' ;atom(22,4)='3d'       
       atom(23,1)='V'  ;atom(23,2)='4s' ;atom(23,3)='4p' ;atom(23,4)='3d'     
       atom(24,1)='Cr' ;atom(24,2)='4s' ;atom(24,3)='4p' ;atom(24,4)='3d'       
       atom(25,1)='Mn' ;atom(25,2)='4s' ;atom(25,3)='4p' ;atom(25,4)='3d'       
       atom(26,1)='Fe' ;atom(26,2)='4s' ;atom(26,3)='4p' ;atom(26,4)='3d'             
       atom(27,1)='Co' ;atom(27,2)='4s' ;atom(27,3)='4p' ;atom(27,4)='3d'                         
       atom(28,1)='Ni' ;atom(28,2)='4s' ;atom(28,3)='4p' ;atom(28,4)='3d'                        
       atom(29,1)='Cu' ;atom(29,2)='4s' ;atom(29,3)='4p' ;atom(29,4)='3d'                        
       atom(30,1)='Zn' ;atom(30,2)='4s' ;atom(30,3)='4p' ;atom(30,4)='3d' 
       atom(31,1)='Ga' ;atom(31,2)='4s' ;atom(31,3)='4p' ;atom(31,4)='3d'
       atom(32,1)='Ge' ;atom(32,2)='4s' ;atom(32,3)='4p' ;atom(32,4)='3d'
       atom(33,1)='As' ;atom(33,2)='4s' ;atom(33,3)='4p' ;atom(33,4)='3d' 
       atom(34,1)='Se' ;atom(34,2)='4s' ;atom(34,3)='4p' ;atom(34,4)='3d'
       atom(35,1)='Br' ;atom(35,2)='4s' ;atom(35,3)='4p' ;atom(35,4)='3d' 
       atom(36,1)='Kr' ;atom(36,2)='4s' ;atom(36,3)='4p' ;atom(36,4)='3d'
      
       atom(37,1)='Rb' ;atom(37,2)='5s' ;atom(37,3)='4p'
       atom(38,1)='Sr' ;atom(38,2)='5s'                  ;atom(38,4)='4d'
       atom(39,1)='Y'  ;atom(39,2)='5s' ;atom(39,3)='5p' ;atom(39,4)='4d'    
       atom(40,1)='Zr' ;atom(40,2)='5s' ;atom(40,3)='5p' ;atom(40,4)='4d'  
       atom(41,1)='Nb' ;atom(41,2)='5s' ;atom(41,3)='5p' ;atom(41,4)='4d'    
       atom(42,1)='Mo' ;atom(42,2)='5s' ;atom(42,3)='5p' ;atom(42,4)='4d'   
       atom(43,1)='Tc' ;atom(43,2)='5s' ;atom(43,3)='5p' ;atom(43,4)='4d'  
       atom(44,1)='Ru' ;atom(44,2)='5s' ;atom(44,3)='5p' ;atom(44,4)='4d' 
       atom(45,1)='Rh' ;atom(45,2)='5s' ;atom(45,3)='5p' ;atom(45,4)='4d'  
       atom(46,1)='Pd' ;atom(46,2)='5s' ;atom(46,3)='5p' ;atom(46,4)='4d'  
       atom(47,1)='Ag' ;atom(47,2)='5s' ;atom(47,3)='5p' ;atom(47,4)='4d'  
       atom(48,1)='Cd' ;atom(48,2)='5s' ;atom(48,3)='5p' ;atom(48,4)='4d'    
       atom(49,1)='In' ;atom(49,2)='5s' ;atom(49,3)='5p' ;atom(49,4)='4d'  
       atom(50,1)='Sn' ;atom(50,2)='5s' ;atom(50,3)='5p' ;atom(50,4)='4d'  
       atom(51,1)='Sb' ;atom(51,2)='5s' ;atom(51,3)='5p' ;atom(51,4)='4d'  
       atom(52,1)='Te' ;atom(52,2)='5s' ;atom(52,3)='5p' ;atom(52,4)='4d'  
       atom(53,1)='I'  ;atom(53,2)='5s' ;atom(53,3)='5p' ;atom(53,4)='4d' 
       atom(54,1)='Xe' ;atom(54,2)='5s' ;atom(54,3)='5p' ;atom(54,4)='4d'  

       atom(55,1)='Cs' ;atom(55,2)='6s' ;atom(55,3)='5p'
       atom(56,1)='Ba' ;atom(56,2)='6s' ;atom(56,3)='5p' ;atom(56,4)='5d'  
       atom(57,1)='La' ;atom(57,2)='6s' ;atom(57,3)='6p' ;atom(57,4)='5d' 
       atom(58,1)='Ce' ;atom(58,2)='6s' ;atom(58,3)='6p'                  ;atom(58,5)='4f' 
       atom(59,1)='Pr' ;atom(59,2)='6s' ;atom(59,3)='6p' ;atom(59,4)='5d' ;atom(59,5)='4f'
       atom(60,1)='Nd' ;atom(60,2)='6s' ;atom(60,3)='6p' ;atom(60,4)='5d' ;atom(60,5)='4f' 

       atom(61,1)='Pm' ;atom(61,2)='6s' ;atom(61,3)='6p' ;atom(61,4)='5d' ;atom(61,5)='4f'                    
       atom(62,1)='Sm' ;atom(62,2)='6s' ;atom(62,3)='6p' ;atom(62,4)='5d' ;atom(62,5)='4f'      
       atom(63,1)='Eu' ;atom(63,2)='6s' ;atom(63,3)='6p' ;atom(63,4)='5d' ;atom(63,5)='4f'                     
       atom(64,1)='Gd' ;atom(64,2)='6s' ;atom(64,3)='6p' ;atom(64,4)='5d' ;atom(64,5)='4f'       
       atom(65,1)='Tb' ;atom(65,2)='6s' ;atom(65,3)='6p' ;atom(65,4)='5d' ;atom(65,5)='4f'       
       atom(66,1)='Dy' ;atom(66,2)='6s' ;atom(66,3)='6p' ;atom(66,4)='5d' ;atom(66,5)='4f'       
       atom(67,1)='Ho' ;atom(67,2)='6s' ;atom(67,3)='6p' ;atom(67,4)='5d' ;atom(67,5)='4f'                        
       atom(68,1)='Er' ;atom(68,2)='6s' ;atom(68,3)='6p' ;atom(68,4)='5d' ;atom(68,5)='4f'                        
       atom(69,1)='Tm' ;atom(69,2)='6s' ;atom(69,3)='6p' ;atom(69,4)='5d' ;atom(69,5)='4f'                       
       atom(70,1)='Yb' ;atom(70,2)='6s' ;atom(70,3)='6p' ;atom(70,4)='5d' ;atom(70,5)='4f'       
       atom(71,1)='Lu' ;atom(71,2)='6s' ;atom(71,3)='6p' ;atom(71,4)='5d' ;atom(71,5)='4f'         
       atom(72,1)='Hf' ;atom(72,2)='6s' ;atom(72,3)='6p' ;atom(72,4)='5d' ;atom(72,5)='4f'         
       atom(73,1)='Ta' ;atom(73,2)='6s' ;atom(73,3)='6p' ;atom(73,4)='5d' ;atom(73,5)='4f'       
       atom(74,1)='W'  ;atom(74,2)='6s' ;atom(74,3)='6p' ;atom(74,4)='5d' ;atom(74,5)='4f'       
       atom(75,1)='Re' ;atom(75,2)='6s' ;atom(75,3)='6p' ;atom(75,4)='5d' ;atom(75,5)='4f'        
       atom(76,1)='Os' ;atom(76,2)='6s' ;atom(76,3)='6p' ;atom(76,4)='5d' ;atom(76,5)='4f'       
       atom(77,1)='Ir' ;atom(77,2)='6s' ;atom(77,3)='6p' ;atom(77,4)='5d' ;atom(77,5)='4f'       
       atom(78,1)='Pt' ;atom(78,2)='6s' ;atom(78,3)='6p' ;atom(78,4)='5d' ;atom(78,5)='4f'     
       atom(79,1)='Au' ;atom(79,2)='6s' ;atom(79,3)='6p' ;atom(79,4)='5d' ;atom(79,5)='4f'       
       atom(80,1)='Hg' ;atom(80,2)='6s' ;atom(80,3)='6p' ;atom(80,4)='5d' ;atom(80,5)='4f'        
       atom(81,1)='Tl' ;atom(81,2)='6s' ;atom(81,3)='6p' ;atom(81,4)='5d' ;atom(81,5)='4f'      
       atom(82,1)='Pb' ;atom(82,2)='6s' ;atom(82,3)='6p' ;atom(82,4)='5d' ;atom(82,5)='4f'  
       atom(83,1)='Bi' ;atom(83,2)='6s' ;atom(83,3)='6p' ;atom(83,4)='5d' ;atom(83,5)='4f'  
       atom(84,1)='Po' ;atom(84,2)='6s' ;atom(84,3)='6p' ;atom(84,4)='5d' ;atom(84,5)='4f'  
       atom(85,1)='At' ;atom(85,2)='6s' ;atom(85,3)='6p' ;atom(85,4)='5d' ;atom(85,5)='4f'  
       atom(86,1)='Rn' ;atom(86,2)='6s' ;atom(86,3)='6p' ;atom(86,4)='5d' ;atom(86,5)='4f'  

       atom(87,1)='Fr' ;atom(87,2) ='7s'
       atom(88,1)='Ra' ;atom(88,2) ='7s'
       atom(89,1)='Ac' ;atom(89,2) ='7s'
       atom(90,1)='Th' ;atom(90,2) ='7s'
       atom(91,1)='Pa' ;atom(91,2) ='7s'                 ;atom(91,4)='6d' ;atom(91,5) ='5f'  
       atom(92,1)='U'  ;atom(92,2) ='7s'                 ;atom(92,4)='6d' ;atom(92,5) ='5f'  
       atom(93,1)='Np' ;atom(93,2) ='7s'                 ;atom(93,4)='6d' ;atom(93,5) ='5f'  
       atom(94,1)='Pu' ;atom(94,2) ='7s'                                  ;atom(94,5) ='5f'
       atom(95,1)='Am' ;atom(95,2) ='7s'                                  ;atom(95,5) ='5f'
       atom(96,1)='Cm' ;atom(96,2) ='7s'                 ;atom(96,4)='6d' ;atom(96,5) ='5f'  
       atom(97,1)='Bk' ;atom(97,2) ='7s'                                  ;atom(97,5) ='5f'
       atom(98,1)='Cf' ;atom(98,2) ='7s'                                  ;atom(98,5) ='5f'
       atom(99,1)='Es' ;atom(99,2) ='7s'                                  ;atom(99,5) ='5f'
       atom(100,1)='Fm';atom(100,2)='7s'                                  ;atom(100,5)='5f'
       atom(101,1)='Md';atom(101,2)='7s'                                  ;atom(101,5)='5f'
       atom(102,1)='No';atom(102,2)='7s'                                  ;atom(102,5)='5f'
       atom(103,1)='Lr';atom(103,2)='7s'                 ;atom(103,4)='6d';atom(103,5)='5f'  
       atom(104,1)='Rf';atom(104,2)='7s'                 ;atom(104,4)='6d';atom(104,5)='5f'  
       atom(105,1)='Db';atom(105,2)='7s'                 ;atom(105,4)='6d';atom(105,5)='5f'  
       atom(106,1)='Sg';atom(106,2)='7s'                 ;atom(106,4)='6d';atom(106,5)='5f'  
       atom(107,1)='Bh';atom(107,2)='7s'                 ;atom(107,4)='6d';atom(107,5)='5f'  
       atom(108,1)='Hs';atom(108,2)='7s'                 ;atom(108,4)='6d';atom(108,5)='5f'  
       atom(109,1)='Mt';atom(109,2)='7s'                 ;atom(109,4)='6d';atom(109,5)='5f'  
       atom(110,1)='Ds';atom(110,2)='7s'                 ;atom(110,4)='6d';atom(110,5)='5f'  
       atom(111,1)='Rg';atom(111,2)='7s'                 ;atom(111,4)='6d';atom(111,5)='5f'  


!Default config __________________________

  ALLOCATE(composition(1:nat,1:5))
  ALLOCATE(reconf(1:nat,1:5))
       reconf = '  '
  
  DO i=1,nat
   DO n=1,111
    IF (aname(i)==atom(n,1))THEN
      DO j=1,5
         composition(i,j)=atom(n,j)
      END DO
    END IF 
   END DO
  END DO
     PRINT *,'Valence orbitals according to periodic table data:'
  DO i=1,nat
     PRINT *,(composition(i,j),j=1,5)
  END DO

 !! Reconfig by user ____________________


  PRINT  *, 'Reconfigure valence orbitals manually (y/N)'
  READ(*,'(a)') answer

  IF (answer.EQ.'Y'.OR.answer.EQ.'y')THEN
         
  !! Write and crate template file: "valence.reconfig" !
 
  WRITE(24,'(A)')'# Configuration of valence orbitals considered in XPS'
  WRITE(24,'(A)')'# Edit this file add/change/remove orbitals. Obey the format!'
  WRITE(24,'(A)')'# Note: You can include only one "s","p","d" or "f" orbital for each atom (eg. 3s or 4s, but NEVER both of them)'
  DO i=1,nat
     counter = 1
     DO n=2,5
     READ(composition(i,n),50) temp1,temp2
     50 FORMAT (A1,A1)
     IF(temp2.EQ.'s'.OR.temp2.EQ.'p'.OR.temp2.EQ.'d'.OR.temp2.EQ.'f')THEN
        counter = counter+1
     END IF   
     END DO      
     WRITE(24,100)'ATOM',i,':',(composition(i,j),j=1,counter)
     100 FORMAT(A4,I5,A1, 5 A2)
          
  END DO
  !END IF
  CLOSE (24)

  CALL system ("$EDITOR valence.reconfig")
  
!!Read data from "valence.reconfig" file 
  
  OPEN (UNIT=114 , FILE='valence.reconfig')
  READ(114,*)
  READ(114,*)
  READ(114,*)
  DO i=1,nat
     read(114,'(a)',iostat=ios) string
     IF (ios.ne.0) goto 120
     do j=1,69
     READ(string(j:j),'(a1)') answer
     if(answer.eq.':') then
       read(string(j+1:j+11),'(5a2)') (reconf(i,lt),lt=1,5)
       goto 121
     endif
     enddo
 121 continue
  END DO 
  CLOSE (114)

  DO i=1,nat
        PRINT *,aname(i),reconf(i,1)
     IF(aname(i).NE.reconf(i,1))THEN
        PRINT *, 'Wrong atom name!'
        ios = -1
        GOTO 120  
     END IF
  END DO



  PRINT*,'____________________________________________________'
  WRITE(*,'(A)')'Default config    Your input config' 
    DO i=1,nat
       WRITE(*,130)(composition(i,j),j=1,5),(reconf(i,j),j=1,5)
       130 FORMAT(A2,A2,A2,A2,A2,8x,A2,A2,A2,A2,A2)
       DO j=1,5
           composition(i,j)=reconf(i,j)
       END DO
    END DO

  120 IF(ios.LT.0)THEN
      PRINT *, 'Incorrect input, continuing with the default valence configuration!'  
      END IF
 END IF
   

END SUBROUTINE periodic_table

SUBROUTINE estimate_cross (database,scheme)

     use struct      ,only : mult, nat, ndif, zz, aname 
     use find_orbital,only : cross_section,multi_cros,mult_spec,stat
     use main        ,only : excitation_energy,Number_of_case_DOS,output_counter,output_dos,trow,output_names,atomic_number,basename,crfname
     use p_table,only      : atom

 IMPLICIT NONE
 INTEGER, INTENT (IN)                 :: database,scheme
 CHARACTER*2                          :: atomname,poname,temp,temp3
 CHARACTER(len=256)                   :: buffer
 INTEGER                              :: ios,io,ls1,temp2
 INTEGER                              :: I,n,counter,temp4,no,l
 INTEGER                              :: orbital_split
 REAL *16                             :: sum_all
 Logical                              :: data_exist = .false.
 INTEGER                              :: atomic_num
 INTEGER ,DIMENSION(2000)             :: atomic
 
 CHARACTER*4                          :: TEMPO,temp6
 CHARACTER*1                          :: answer
 INTEGER                              :: str,stp,J,k

 REAL*16 , DIMENSION  (2000)          :: atomic_z   
 REAL*16 , DIMENSION  (2000)          :: calc_cross 

 INTEGER                              :: temp5,temp7,es_counter
 CHARACTER*10,dimension (1:1000)      :: es_names1
 CHARACTER*2 ,dimension (1:1000)      :: es_names2
 REAL*8 , DIMENSION (1:118,1:30)      :: es_cross
 CHARACTER*180                        :: esfname
 character*2                          :: cnumber
!################################################################
!! Reading Data from Database2 in order to craete Template file for estimation!!
  ios        = 0
  es_counter = 0
  es_cross   = 0.0
  WRITE(23,'(A)')'# Template file _ generated for the estimation of atomic cross sections'
  WRITE(23,'(A)')'# Edit the file in order to include or exclude orbital name'
  WRITE(23,'(A)')'# Follow the format '

  DO J=1,output_counter
     IF (cross_section(J).EQ.0.0) THEN
         READ(output_names(J),701) temp6
         read(OUTPUT_NAMES(j)(5:9),*) temp5
         701 FORMAT (A4) 
       
             DO I=1,output_counter
              READ(output_names(I),701) temp6
              READ(output_names(I)(5:9),*) temp7
                IF(temp5.EQ.temp7)THEN
                    es_counter = es_counter+1
                    READ(output_names(I),702)  es_names2(es_counter)
                    702 FORMAT (2x,A2)
                    es_names1(es_counter) = output_names(I)
                
                      DO k=1,es_counter-1
                        IF(es_names2(es_counter).EQ.es_names2(K))THEN
                           es_counter = es_counter-1
                        END IF    
                      END DO

                END IF
             END DO

      END IF
 END DO 
!#####################################
! Write data into template file '  
 DO I=1,es_counter 
    WRITE(23,703) es_names2(I) 
    703 FORMAT (A2)
 END DO
 CLOSE (23)
!#####################################
! Open template file for user '
  CALL system ("$EDITOR pes.template")
!#####################################  
! Read input data from template file ' 

 OPEN (UNIT=114 , FILE='pes.template')
       READ(114,*)
       READ(114,*)
       READ(114,*)
       es_counter = 0
 DO WHILE (ios .EQ. 0) 
     es_counter=es_counter+1
     READ(114, '(A)', IOSTAT=ios) es_names1(es_counter)
       IF (ios.NE.0) THEN
         es_counter=es_counter-1
       END IF  
 END DO 
 CLOSE (114)


 IF (es_counter.EQ.0) THEN
    GOTO 1000
 END IF


!#####################################
!Colculate requsted partial orbital cross section from DATABASE1

IF (database .EQ. 1 ) THEN 
 
DO J = 1,es_counter
   READ(4,*)
   ios = 0
   str = 1
   stp = 2 
   READ(es_names1(J),100)poname
   100 FORMAT (A2)
   counter = 0
   DO WHILE (ios == 0)
      READ(4,'(A)', IOSTAT=ios, END =50) buffer
      50 ls1 = len_trim(buffer)
         IF(ls1.NE.0) THEN 
            READ (buffer,*,END=60) temp
      60 IF (temp.EQ.poname) THEN
             BACKSPACE (4)
             BACKSPACE (4)
             READ(4,200)temp3,atomic_num
             200 FORMAT(A2,3x,I3)         
               counter          = counter +1 
               no               = output_counter+1+counter
               TEMPO            = temp3 // poname
               atomic (no)      = atomic_num
               output_names(no) = TEMPO
                  IF (atomic (no).EQ.atomic (no-1)) THEN
                    counter = counter -1
                  END IF
                  READ (4,*)
                  READ (4,*) 
          END IF
          END IF
         temp='  '
   END DO
   REWIND (4)
   ios = 0

     DO i=1,counter
         no =       output_counter+1+i
         CALL       Read_database2 (database,no,scheme)
         atomic_z   (i) = atomic (no)
         calc_cross (i) = cross_section(no)
        DO k=1,120
         IF (K.EQ.atomic_z(i)) THEN
         es_cross(k,J) =  calc_cross (i)
         END IF
        END DO    
      END DO
          IF (counter.EQ.0) THEN
          PRINT *, 'Database contains no data on:',es_names1(J)                    
          END IF 
END DO


!#####################################
!Colculate requsted partial orbital cross section from DATABASE2

 ELSE IF (database.EQ.2 ) THEN
         
    DO I=1,es_counter
    DO J=1,111
       no = output_counter+1+J
       output_names(no)=  atom(J,1)// es_names1(I)
       CALL Read_database (database,no,basename,2)
            es_cross(J,I)= cross_section(no)
    END DO
    END DO 
END IF
  
!#####################################
!Cross sections output writing
 

  temp7=ABS((es_counter+1)/7) 
  str=1
  stp=6

 DO k=0,temp7
    IF (stp.GT.es_counter) stp=es_counter
    IF(k.LE.1) THEN
       WRITE(cnumber,'(i1)') k+1
       esfname=trim(crfname)//trim(cnumber)
       PRINT*,'generating ',trim(esfname)
       print*, 'examine this file to extrapolate the missing cross section'
       print*, 'enter this estimated cross section later on.'
       OPEN(51+k,FILE=esfname)
    END IF
    IF(k.GT.1) THEN
        IF(k.LT.9) write(cnumber,'(i1)') k+1
        IF(k.GE.9) write(cnumber,'(i2)') k+1
        esfname=trim(crfname)//trim(cnumber)
        PRINT*,'generating ',trim(esfname)
       OPEN(51+k,FILE=esfname)
    END IF
           
    WRITE((51+k),450) '#','Tabulated Cross sections for'
    WRITE((51+k),460) '#','Excitation enegy:',excitation_energy,' eV'
    WRITE((51+k),500)'#',' Atom   nuclear charge',(es_names1(j),j=str,stp)
    
    DO I=1,111
       sum_all = 0.0
       DO l=1,es_counter
          sum_all= sum_all+  es_cross(I,l)
       END DO
           IF(sum_all .GT. 0.0 ) THEN
              WRITE (51+k,600) atom(I,1),I,(es_cross(I,j),j=str,stp)
            
           END IF
    END DO
    str=str+6
    stp=stp+6
    CLOSE (UNIT=51+k)  
    buffer="$EDITOR " // esfname
    CALL system (buffer)
 END DO 

 450 FORMAT (A,2x,A24)  
 460 FORMAT (A,2x,A17,f16.5,A3)
 500 FORMAT (A,A22,7(A16,2x))
 600 FORMAT (2x,a4,4x,I5,5x,7(e16.8,2x))

 1000 END SUBROUTINE


_______________________________________________
Wien mailing list
Wien@zeus.theochem.tuwien.ac.at
http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
SEARCH the MAILING-LIST at:  
http://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/index.html

Reply via email to