|
1 | 1 | program optimwidths |
| 2 | + use m_read |
| 3 | + use m_misc |
2 | 4 | use m_func |
3 | 5 | use m_optimizers |
4 | 6 | implicit none |
@@ -35,7 +37,7 @@ program optimwidths |
35 | 37 | atomTypes_a = 0 |
36 | 38 | call readFreqOut(nrAtoms, cartDOF, normDOF, W, M, U, & |
37 | 39 | masses, prog, fileName, atomTypes_a,& |
38 | | - ntype, constrained) |
| 40 | + ntype) |
39 | 41 | call convert2nm(nrAtoms, cartDOF, normDOF, U, M, masses) |
40 | 42 | call pfunc%mf_init(nrAtoms, cartDOF, normDOF, W, M, U, & |
41 | 43 | atomTypes_a, constrained) |
@@ -67,294 +69,3 @@ program optimwidths |
67 | 69 | deallocate(W, M, A, U, widths, masses, atomTypes_a, & |
68 | 70 | & pfunc, x0) |
69 | 71 | end program optimwidths |
70 | | - |
71 | | -subroutine readInput(nrAtoms, maxiter, xtol, ftol, & |
72 | | - constrained, prog, fName) |
73 | | - integer, intent(inout) :: maxiter, nrAtoms |
74 | | - real(kind=8), intent(inout) :: xtol, ftol |
75 | | - logical, intent(inout) :: constrained |
76 | | - character(len=128), intent(inout) :: prog, fName |
77 | | - namelist / optimparams / nrAtoms, maxiter, xtol, & |
78 | | - ftol, constrained, prog,& |
79 | | - fName |
80 | | - |
81 | | - open(10, file='input.in') |
82 | | - read(unit=10, nml=optimparams) |
83 | | - write(*, '(A)') 'Read user input:' |
84 | | - write(unit=*, nml=optimparams) |
85 | | - close(10) |
86 | | -end subroutine readInput |
87 | | - |
88 | | -subroutine readFreqOut(nrAtoms, cDOF, nDOF, W, M, U, masses, & |
89 | | - prog, fName, atomTypes_a, ntype) |
90 | | - use m_func |
91 | | - implicit none |
92 | | - integer, intent(in) :: nrAtoms, cDOF, nDOF |
93 | | - integer, intent(inout) :: ntype |
94 | | - real(kind=8), intent(inout) :: W(nDOF,nDOF), M(nDOF,nDOF), & |
95 | | - U(cDOF,nDOF), masses(nrAtoms), & |
96 | | - atomTypes_a(nrAtoms) |
97 | | - character(len=128), intent(in) :: prog, fName |
98 | | - character(len=256) :: line |
99 | | - real(kind=8), allocatable :: atomCoords(:,:) |
100 | | - integer, allocatable :: nrMasses(:) |
101 | | - character(len=128), allocatable :: atomNames(:) |
102 | | - character(len=20), allocatable :: tempCLine(:) |
103 | | - real(kind=8), allocatable :: tempRLine(:) |
104 | | - real(kind=8), parameter :: freqConv = 4.5563352529119736d-6 |
105 | | - real(kind=8), parameter :: proton = 1822.888515d0 |
106 | | - logical :: rMasses, rAtomNames, rAtomCoords, rNormalModes |
107 | | - integer :: ioerr, ind, lineNr, i, j, ismass, ifreq, irmass, & |
108 | | - inmode, jnmode, itmp, jtmp, tmpNrAtoms, nMass, & |
109 | | - imass |
110 | | - |
111 | | - |
112 | | - allocate(atomCoords(nrAtoms,3), atomNames(nrAtoms)) |
113 | | - nMass = 1 |
114 | | - tmpNrAtoms = nrAtoms |
115 | | - if (tmpNrAtoms > 10) then |
116 | | - do |
117 | | - tmpNrAtoms = tmpNrAtoms - 10 |
118 | | - nMass = nMass + 1 |
119 | | - if (tmpNrAtoms < 10) exit |
120 | | - enddo |
121 | | - endif |
122 | | - allocate(nrMasses(nMass)) |
123 | | - do i=1,nMass |
124 | | - if (i < nMass) nrMasses(i) = 10 |
125 | | - if (i == nMass) nrMasses(i) = tmpNrAtoms |
126 | | - enddo |
127 | | - rAtomNames = .false. |
128 | | - rAtomCoords = .false. |
129 | | - rNormalModes = .false. |
130 | | - lineNr = 0 |
131 | | - ismass = 1 |
132 | | - imass = 1 |
133 | | - ifreq = 1 |
134 | | - irmass = 1 |
135 | | - jnmode = 0 |
136 | | - inmode = 1 |
137 | | - W = 0.d0 |
138 | | - M = 0.d0 |
139 | | - U = 0.d0 |
140 | | - select case(prog) |
141 | | - case('gaussian') |
142 | | - open(unit=11, file=fName, status='old', action='read', & |
143 | | - iostat=ioerr) |
144 | | - if (ioerr == 0) then |
145 | | - do |
146 | | - read(11, '(A)', iostat=ioerr) line |
147 | | - if (ioerr /= 0) exit |
148 | | - ind = index(trim(line), 'Symbolic Z-matrix') |
149 | | - if (ind /= 0) then |
150 | | - rAtomNames = .true. |
151 | | - lineNr = lineNr + 1 |
152 | | - cycle |
153 | | - endif |
154 | | - |
155 | | - if (rAtomNames .and. (lineNr > 1)) then |
156 | | - read(line, '(x,a2,19x,3f9.5)') atomNames(lineNr-1), & |
157 | | - atomCoords(lineNr-1,:) |
158 | | - atomCoords = 0.d0 |
159 | | - lineNr = lineNr + 1 |
160 | | - if (lineNr > (nrAtoms + 1)) then |
161 | | - lineNr = 0 |
162 | | - rAtomNames = .false. |
163 | | - call fillAtomTypes(nrAtoms,ntype,atomTypes_a,& |
164 | | - atomNames) |
165 | | - endif |
166 | | - elseif (rAtomNames .and. (lineNr == 1)) then |
167 | | - lineNr = lineNr + 1 |
168 | | - endif |
169 | | - |
170 | | - ind = index(trim(line), 'AtmWgt=') |
171 | | - if (ind /= 0) then |
172 | | - if (ismass < nMass) allocate(tempCLine(1),tempRLine(10)) |
173 | | - if (ismass == nMass) allocate(tempCLine(1), & |
174 | | - tempRLine(nrMasses(ismass))) |
175 | | - if (ismass > nMass) cycle |
176 | | - read(line,*) tempCLine(:), tempRLine(:) |
177 | | - if (ismass < nMass) then |
178 | | - do i=1,10 |
179 | | - masses(imass) = tempRLine(i)*proton |
180 | | - !write(*,*) masses(imass) |
181 | | - imass = imass + 1 |
182 | | - enddo |
183 | | - elseif (ismass == nMass) then |
184 | | - do i=1,nrMasses(ismass) |
185 | | - masses(imass) = tempRLine(i)*proton |
186 | | - !write(*,*) masses(imass) |
187 | | - imass = imass + 1 |
188 | | - enddo |
189 | | - endif |
190 | | - ismass = ismass + 1 |
191 | | - deallocate(tempCLine, tempRLine) |
192 | | - endif |
193 | | - |
194 | | - ind = index(trim(line), 'Standard orientation:') |
195 | | - if (ind /= 0) then |
196 | | - rAtomCoords = .true. |
197 | | - lineNr = lineNr + 1 |
198 | | - cycle |
199 | | - endif |
200 | | - |
201 | | - if (rAtomCoords .and. (lineNr > 4)) then |
202 | | - allocate(tempRLine(6)) |
203 | | - read(line, *) tempRLine |
204 | | - atomCoords(lineNr-4,:) = tempRLine(4:) |
205 | | - deallocate(tempRLine) |
206 | | - lineNr = lineNr + 1 |
207 | | - if (lineNr > (nrAtoms + 4)) then |
208 | | - lineNr = 0 |
209 | | - rAtomCoords = .false. |
210 | | - endif |
211 | | - elseif (rAtomCoords .and. (lineNr <= 4)) then |
212 | | - lineNr = lineNr + 1 |
213 | | - endif |
214 | | - |
215 | | - ind = index(trim(line), 'Frequencies') |
216 | | - if (ind /= 0) then |
217 | | - rNormalModes = .true. |
218 | | - lineNr = lineNr + 1 |
219 | | - allocate(tempCLine(2),tempRLine(3)) |
220 | | - read(line,*) tempCLine(:), tempRLine(:) |
221 | | - do i=1,3 |
222 | | - W(ifreq,ifreq) = tempRLine(i)*freqConv |
223 | | - ifreq = ifreq + 1 |
224 | | - enddo |
225 | | - deallocate(tempCLine,tempRLine) |
226 | | - cycle |
227 | | - endif |
228 | | - |
229 | | - if (rNormalModes .and. (lineNr == 1)) then |
230 | | - allocate(tempCLine(3),tempRLine(3)) |
231 | | - read(line,*) tempCLine(:), tempRLine(:) |
232 | | - do i=1,3 |
233 | | - M(irmass,irmass) = tempRLine(i)*proton |
234 | | - irmass = irmass + 1 |
235 | | - enddo |
236 | | - deallocate(tempCLine,tempRLine) |
237 | | - lineNr = lineNr + 1 |
238 | | - elseif (rNormalModes .and. (lineNr <= 4)) then |
239 | | - lineNr = lineNr + 1 |
240 | | - elseif (rNormalModes .and. (lineNr > 4)) then |
241 | | - lineNr = lineNr + 1 |
242 | | - allocate(tempRLine(2)) |
243 | | - read(line,*) tempRLine(:), & |
244 | | - U(inmode:inmode+2,jnmode+1), & |
245 | | - U(inmode:inmode+2,jnmode+2), & |
246 | | - U(inmode:inmode+2,jnmode+3) |
247 | | - inmode = inmode + 3 |
248 | | - deallocate(tempRLine) |
249 | | - if (lineNr > (nrAtoms + 4)) then |
250 | | - jnmode = jnmode + 3 |
251 | | - inmode = 1 |
252 | | - lineNr = 0 |
253 | | - rNormalModes = .false. |
254 | | - endif |
255 | | - endif |
256 | | - enddo |
257 | | - endif |
258 | | - close(11) |
259 | | - endselect |
260 | | - deallocate(atomCoords, atomNames) |
261 | | -end subroutine readFreqOut |
262 | | - |
263 | | -subroutine convert2nm(nrAtoms, cDOF, nDOF, U, M, masses) |
264 | | - integer, intent(in) :: nrAtoms, cDOF, nDOF |
265 | | - real(kind=8), intent(in) :: M(nDOF,nDOF), masses(nrAtoms) |
266 | | - real(kind=8), intent(inout) :: U(cDOF,nDOF) |
267 | | - integer :: i, j, iAtom |
268 | | - |
269 | | - do i = 1,cDOF |
270 | | - iAtom = (i-1)/3 + 1 |
271 | | - do j = 1,nDOF |
272 | | - U(i,j) = (U(i,j)*dsqrt(masses(iAtom))/dsqrt(M(j,j))) |
273 | | - enddo |
274 | | - enddo |
275 | | - |
276 | | -end subroutine convert2nm |
277 | | - |
278 | | -subroutine fillAtomTypes(nrAtoms, ind, atomTypes_a, atomNames) |
279 | | - use m_func |
280 | | - implicit none |
281 | | - integer, intent(in) :: nrAtoms |
282 | | - integer, intent(inout) :: ind |
283 | | - character(len=128), intent(in) :: atomNames(nrAtoms) |
284 | | - integer, intent(inout) :: atomTypes_a(nrAtoms) |
285 | | - type(mf_llist), pointer :: head => NULL(), tail => NULL(), & |
286 | | - ptr => NULL(), ptr1 => NULL(), & |
287 | | - next => NULL(), current => NULL() |
288 | | - integer :: i |
289 | | - logical :: inllist |
290 | | - |
291 | | - allocate(ptr) |
292 | | - ind = ind + 1 |
293 | | - do i=1,nrAtoms |
294 | | - inllist = .false. |
295 | | - ptr%atomName = atomNames(i) |
296 | | - ptr%atomIndex = ind |
297 | | - if (not(associated(head))) then |
298 | | - allocate(head) |
299 | | - tail => head |
300 | | - nullify(ptr%mf_next) |
301 | | - tail%atomName = ptr%atomName |
302 | | - tail%atomIndex = ptr%atomIndex |
303 | | - ind = ind + 1 |
304 | | - elseif (ptr%atomName /= tail%atomName) then |
305 | | - ptr1 => head |
306 | | - do |
307 | | - if (ptr%atomName == ptr1%atomName) then |
308 | | - inllist = .true. |
309 | | - exit |
310 | | - endif |
311 | | - ptr1 => ptr1%mf_next |
312 | | - if (not(associated(ptr1))) exit |
313 | | - enddo |
314 | | - if (not(inllist)) then |
315 | | - allocate(tail%mf_next) |
316 | | - tail => tail%mf_next |
317 | | - nullify(tail%mf_next) |
318 | | - tail%atomName = ptr%atomName |
319 | | - tail%atomIndex = ptr%atomIndex |
320 | | - ind = ind + 1 |
321 | | - endif |
322 | | - endif |
323 | | - enddo |
324 | | - deallocate(ptr) |
325 | | - |
326 | | - ind = ind - 1 |
327 | | - write(*,*) "There are ", ind, " distinct atom types:" |
328 | | - write(*,'(A,10x,A)') "Name", "Index" |
329 | | - ptr => head |
330 | | - do |
331 | | - write(*,*) ptr%atomName, ptr%atomIndex |
332 | | - ptr => ptr%mf_next |
333 | | - if(not(associated(ptr))) exit |
334 | | - enddo |
335 | | - |
336 | | - do i=1,nrAtoms |
337 | | - ptr => head |
338 | | - do |
339 | | - if (ptr%atomName == atomNames(i)) then |
340 | | - atomTypes_a(i) = ptr%atomIndex |
341 | | - endif |
342 | | - ptr => ptr%mf_next |
343 | | - if(not(associated(ptr))) exit |
344 | | - enddo |
345 | | - enddo |
346 | | - |
347 | | - current => head |
348 | | - next => current%mf_next |
349 | | - do |
350 | | - deallocate(current) |
351 | | - nullify(current) |
352 | | - if(not(associated(next))) exit |
353 | | - current => next |
354 | | - next => current%mf_next |
355 | | - enddo |
356 | | - |
357 | | - nullify(head) |
358 | | - nullify(tail) |
359 | | - |
360 | | -end subroutine fillAtomTypes |
0 commit comments