Statistics
| Revision:

root / src / ReadInput_mopac.f90 @ 10

History | View | Annotate | Download (4.8 kB)

1
 SUBROUTINE ReadInput_Mopac
2

    
3
! This routine reads an input template for MOPAC
4

    
5
  use VarTypes
6
  use Path_module
7
  use Io_module
8

    
9
  IMPLICIT NONE
10

    
11
  INTERFACE
12
     function valid(string) result (isValid)
13
       CHARACTER(*), intent(in) :: string
14
       logical                  :: isValid
15
     END function VALID
16

    
17
    SUBROUTINE die(routine, msg, file, line, unit)
18

    
19
      Use VarTypes
20
      Use io_module
21

    
22
      implicit none
23
      character(len=*), intent(in)           :: routine, msg
24
      character(len=*), intent(in), optional :: file
25
      integer(KINT), intent(in), optional      :: line, unit
26

    
27
    END SUBROUTINE die
28

    
29

    
30
  END INTERFACE
31

    
32

    
33
  CHARACTER(LCHARS) ::  Line,LineUp
34
  INTEGER(KINT) :: LineL, Idx, NTmp
35
  INTEGER(KINT) :: NatMopac
36
  REAL(KREAL) :: Lat(3,3)
37
  
38
  LOGICAL :: Debug
39

    
40

    
41
  Debug=Valid("readinput").OR.Valid("readinput_mopac")
42

    
43
  if (debug) Call Header("Entering ReadInput_mopac")
44

    
45
! The structure is:
46
! A MOPAC data set normally consists of one line of keywords, two lines of user-defined text, then the coordinates
47
! Then  a blank line or a line of 0.
48
! then the symmetry description.
49
! comment lines start with * and can be anywhere !!!
50

    
51
  ! First, the root
52
  IF (DEBUG) WRITE(*,*) "Reading Mopac input"
53
  ALLOCATE(Mopac_Root)
54
  NULLIFY(Mopac_Root%next)
55
  ALLOCATE(Mopac_Comment)
56
  NULLIFY(Mopac_Comment%next)
57
  ALLOCATE(Mopac_End)
58
  NuLLIFY(Mopac_End%Next)
59
  Current => Mopac_root
60
  CurCom => Mopac_Comment
61
  LineL=1
62
  NTmp=0
63
  DO WHILE (NTmp.LT.3)
64
     READ(IOIN,'(A)') Line
65
     Line=AdjustL(Line)
66
     LineL=len(Trim(Line))
67
     IF (Line(1:1)/="*") THEN
68
        IF (NTmp==0) THEN
69
           LineUp=Line
70
           Call UpCase(LineUp)
71
           Idx=Index(LineUp,'GRADIENTS')
72
           If (Idx==0) Line=TRIM(Line) // " GRADIENTS"
73
           Idx=Index(LineUp,'1SCF')
74
           If (Idx==0) Line=TRIM(Line) // " 1SCF"
75
        END IF
76
        current%Line=TRIM(Line)
77
        ALLOCATE(current%next)
78
        Current => Current%next
79
        Nullify(Current%next)
80
        NTmp=NTmp+1
81
     ELSE
82
        CurCom%Line=TRIM(LINE)
83
        ALLOCATE(CurCom%Next)
84
        CurCom => CurCom%Next
85
        NULLIFY(CurCom%Next)
86
     END IF
87
  END DO
88

    
89
!     Current => Mopac_root
90
!     DO WHILE (ASSOCIATED(Current%next))
91
!        WRITE(*,'(1X,A)') Trim(current%line)
92
!        Current => current%next
93
!     END DO
94

    
95
! Now the geometry... that we just skip :)
96
! PFL 2013 Apr 
97
! We take care that there is no Translation vectors...
98
! We also check that the number of atoms is ok
99
  IF (DEBUG) WRITE(*,*) "Reading Mopac Geometry"
100
  Mopac_EndGeom=""
101
  LineL=1
102
  NatMopac=0
103
  Lat=0.d0
104
  IPer=0
105
  FPBC=.FALSE.
106
  DO WHILE (LineL.NE.0)
107
     READ(IOIN,'(A)',END=989) Line
108
     Line=AdjustL(Line)
109
     LineL=len(Trim(Line))
110
     ! The last line might be either blank or filled with 0
111
     If (LineL>0) THEN
112
        SELECT CASE (Line(1:1))
113
          CASE ("0") 
114
             LineL=0
115
             Mopac_EndGeom=Trim(Line)
116
          CASE("*") 
117
             CurCom%Line=TRIM(LINE)
118
             ALLOCATE(CurCom%Next)
119
             CurCom => CurCom%Next
120
             NULLIFY(CurCom%Next)
121
          CASE DEFAULT
122
             LineUp=Line
123
             Call UpCase(LineUp)
124
             If (LineUp(1:2)=="TV") THEN
125
                FPBC=.TRUE.
126
                IPer=IPer+1
127
                If (Iper>3) THEN
128
                   Call Die("ReadInput Mopac","Iper>3",Unit=IOOUT)
129
                END IF
130
                NTmp=Index(LineUp," ")
131
                LineUp=LineUp(NTmp:)
132
                Read(LineUp,*) Lat(IPer,1:3)
133
             ELSE
134
                NatMopac=NatMopac+1
135
             END IF
136
         END SELECT
137
      END IF
138
     
139
   END DO
140

    
141
!  WRITE(*,*) "NatMopac,Nat:",NAtMopac,Nat
142
  IF (NatMopac/=Nat) Call Die("ReadInput_mopac","Nat read does not mat nat",Unit=IOOUT)
143
  IF (FPBC) THEN
144
     Lat_a(1:3)=Lat(1,1:3)
145
     Lat_b(1:3)=Lat(2,1:3)
146
     Lat_c(1:3)=Lat(3,1:3)
147
     If (IPer>=1) THEN
148
        kaBeg=-1
149
        kaEnd=1
150
     END IF
151
     If (IPer>=2) THEN
152
        kbBeg=-1
153
        kbEnd=1
154
     END IF
155
     If (IPer==3) THEN
156
        kcBeg=-1
157
        kcEnd=1
158
     END IF
159
     If (IPer>3) THEN
160
        Call Die("Readinput_mopac","Found too many Tv lines !",Unit=IOOUT)
161
     END IF
162
  END IF
163

    
164
! If we are here, there might be something else to read: Mopac_end
165

    
166
  ! We now read the last part
167
  IF (DEBUG) WRITE(*,*) "Reading Mopac End"
168
  !     READ(IOIN,'(A)') Line
169
  Current => Mopac_End
170
  LineL=1
171
  DO WHILE (1.EQ.1)
172
     READ(IOIN,'(A)',END=989) Line
173
     Line=AdjustL(Line)
174
     LineL=len(Trim(Line))
175
     IF (Line(1:1)/="*") THEN
176
        current%Line=TRIM(Line)
177
        ALLOCATE(current%next)
178
        Current => Current%next
179
        Nullify(Current%next)
180
        NTmp=NTmp+1
181
     ELSE
182
        CurCom%Line=TRIM(LINE)
183
        ALLOCATE(CurCom%Next)
184
        CurCom => CurCom%Next
185
        NULLIFY(CurCom%Next)
186
        END IF
187
     END DO
188
989  CONTINUE
189

    
190
  if (debug) Call Header("Exiting ReadInput_mopac")
191

    
192
END SUBROUTINE READINPUT_Mopac