← Back to branch summary

~ubuntu-branches/ubuntu/utopic/nwchem/utopic

« back to all changes in this revision

Viewing changes to src/nwdft/rt_tddft/rtutils/rt_tddft_movecs_import.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/nwdft/rt_tddft/rtutils/rt_tddft_movecs_import.F
 
1
C====================================================================
 
2
C
 
3
C     Determine name of initial movecs file.  If one was not specified
 
4
C     in the input deck, it falls back to SCF output.
 
5
C
 
6
      subroutine rt_tddft_movecs_fname (params, rt_movecs_fname)
 
7
 
 
8
      implicit none
 
9
 
 
10
#include "errquit.fh"
 
11
#include "mafdecls.fh"
 
12
#include "stdio.fh"
 
13
#include "global.fh"
 
14
#include "msgids.fh"
 
15
#include "geom.fh"
 
16
#include "util.fh"
 
17
#include "cdft.fh"
 
18
#include "rtdb.fh"
 
19
#include "rt_tddft.fh"
 
20
 
 
21
 
 
22
C     == Inputs ==
 
23
      type(rt_params_t), intent(in) :: params
 
24
 
 
25
      
 
26
C     == Outputs ==
 
27
      character(len=*), intent(out) :: rt_movecs_fname
 
28
 
 
29
 
 
30
C     == Parameters ==
 
31
      character(*),parameter :: pname="rt_tddft_movecs_fname: "
 
32
 
 
33
 
 
34
 
 
35
C
 
36
C     If we found a "vectors" directive while parsing the input we use
 
37
C     that, otherwise fall back on the output vectors "movecs_in" from
 
38
C     the DFT (via common block in cdft.fh header).
 
39
C
 
40
      if (.not. rtdb_cget (params%rtdb, "rt_tddft:init_movecs",
 
41
     $     1, rt_movecs_fname)) then
 
42
 
 
43
         if (movecs_in .eq. "atomic")
 
44
     $        call errquit ('Cannot use atomic guess for rt_tddft'//
 
45
     $        ' initial state--first run an SCF to generate'//
 
46
     $        ' valid movecs', 0, 0)
 
47
 
 
48
         call rt_tddft_print_warning ("Starting movecs not specified"//
 
49
     $        "--trying SCF output: "//char(10)//"   "//
 
50
     $        trim(movecs_out))
 
51
 
 
52
         rt_movecs_fname = movecs_in
 
53
      endif
 
54
 
 
55
      call util_file_name_resolve (rt_movecs_fname, .false.)
 
56
 
 
57
      end subroutine
 
58
 
 
59
 
 
60
 
 
61
C====================================================================
 
62
C
 
63
C     Print imported movecs information to stdout.
 
64
C
 
65
      subroutine rt_tddft_movecs_print_header (params, fname,
 
66
     $     title, basis_name, scftype, nbf, mo_nsets, nmo)
 
67
 
 
68
      implicit none
 
69
 
 
70
#include "errquit.fh"
 
71
#include "stdio.fh"
 
72
#include "global.fh"
 
73
#include "util.fh"
 
74
#include "rt_tddft.fh"
 
75
 
 
76
 
 
77
C     == Inputs ==
 
78
      type(rt_params_t), intent(in) :: params
 
79
      character(len=*), intent(in)  :: fname, title, basis_name, scftype
 
80
      integer, intent(in)           :: nbf, mo_nsets, nmo(mo_nsets)
 
81
 
 
82
 
 
83
C     == Parameters ==
 
84
      character(*),parameter :: pname="rt_tddft_movecs_print_header: "
 
85
 
 
86
 
 
87
C     == Variables ==
 
88
      integer me
 
89
 
 
90
      
 
91
      me = ga_nodeid ()
 
92
 
 
93
      if (me.eq.0) then
 
94
         write (luout, "(5x,a,a)")  "File name              : ",
 
95
     $        trim (fname)
 
96
         write (luout, "(5x,a,a)")  "Job title              : ",
 
97
     $        trim (title)
 
98
         write (luout, "(5x,a,a)")  "Basis set name         : ",
 
99
     $        trim (basis_name)
 
100
         write (luout, "(5x,a,a)")  "SCF type               : ",
 
101
     $        trim (scftype)
 
102
         write (luout, "(5x,a,i0)") "Atomic orbitals        : ",
 
103
     $        nbf
 
104
 
 
105
         if (mo_nsets .eq. 1) then
 
106
            write (luout, "(5x,a,i0)") "Molecular orbitals     : ",
 
107
     $           nmo
 
108
         elseif (mo_nsets .eq. 2) then
 
109
            write (luout, "(5x,a,i0,a,i0)") "Molecular orbitals     : ",
 
110
     $           nmo(1), ", ", nmo(2)
 
111
         endif
 
112
            
 
113
         call util_flush (luout)
 
114
      endif
 
115
      call ga_sync ()
 
116
 
 
117
      end subroutine
 
118
 
 
119
 
 
120
C====================================================================
 
121
C
 
122
C     Print imported movecs eigenvalues and occupations to stdout.
 
123
C
 
124
      subroutine rt_tddft_movecs_print_evals (params, occ, evals)
 
125
 
 
126
      implicit none
 
127
 
 
128
#include "errquit.fh"
 
129
#include "stdio.fh"
 
130
#include "global.fh"
 
131
#include "util.fh"
 
132
#include "rt_tddft.fh"
 
133
 
 
134
 
 
135
C     == Inputs ==
 
136
      type(rt_params_t), intent(in) :: params
 
137
      double precision, intent(in)  :: occ(*), evals(*)
 
138
 
 
139
 
 
140
C     == Parameters ==
 
141
      character(*),parameter :: pname="rt_tddft_movecs_print_evals: "
 
142
 
 
143
 
 
144
C     == Variables ==
 
145
      integer me
 
146
      integer iv
 
147
 
 
148
      
 
149
      me = ga_nodeid ()
 
150
 
 
151
C     
 
152
C     Print eigenvalue list
 
153
C     
 
154
      if (me.eq.0) then
 
155
            
 
156
         write (luout, *) ""
 
157
            
 
158
         call util_print_centered (luout, "  Vector      "//
 
159
     $        "Occupation     Eigenvalue [au]", 0, .true.)
 
160
         
 
161
         do iv = 1, params%ns_ao  
 
162
            if (iv .le. params%ns_mo) then ! vals not removed by canorg
 
163
               write (luout, "(2x,i8,7x,1f5.2,7x,es15.8)")
 
164
     $              iv, occ(iv), evals(iv)
 
165
            else ! lindep vecs
 
166
               write (luout, "(2x,i8,4x,a)")
 
167
     $              iv, "  (linearly dependent; removed)  "
 
168
            endif
 
169
         enddo
 
170
         
 
171
         write (luout, *) ""
 
172
         call util_flush (luout)
 
173
      endif
 
174
      call ga_sync ()
 
175
 
 
176
      end subroutine
 
177
 
 
178
 
 
179
C====================================================================
 
180
C
 
181
C     Convert real-valued ns_ao x ns_ao into complex ns_ao x ns_ao dens
 
182
C     mat in AO basis (used for importing from the ground state code).
 
183
C
 
184
C     P_uv = \sum_i n_i C_ui C_vi^*
 
185
C
 
186
C     where n_i is the ith orbital occupation, and C is the movecs
 
187
C     coefficient matrix.
 
188
C
 
189
C     We first build mask of orbital occupations (n3 = occupation of 3rd
 
190
C     MO).
 
191
C
 
192
C     [ n1   n2  ... nN ]
 
193
C     [ n1   n2      nN ]
 
194
C     [ ...  ...     ...]
 
195
C     [ n1   n2      nN ]
 
196
C
 
197
C     We apply this mask to the movecs matrix (g_movecs_pad) then take
 
198
C     matrix multiplication with the unmasked transpose (note: full
 
199
C     matrices are multiplied, since the mask has zeroed out the
 
200
C     non-occupied orbitals anyways).
 
201
C     
 
202
      subroutine rt_tddft_movecs_zdens(params, nsets, g_zdens_ao)
 
203
 
 
204
      implicit none
 
205
 
 
206
#include "errquit.fh"
 
207
#include "stdio.fh"
 
208
#include "global.fh"
 
209
#include "mafdecls.fh"
 
210
#include "util.fh"
 
211
#include "rt_tddft.fh"
 
212
 
 
213
 
 
214
C     == Inputs ==
 
215
      type(rt_params_t), intent(in) :: params
 
216
C      integer, intent(in)           :: ispin  !target movecs set    : 1 for CS/SO; 1 or 2 for OS
 
217
      integer, intent(in)           :: nsets  !total num movecs sets: 1 for CS/SO; 2 for OS
 
218
 
 
219
 
 
220
C     == Outputs ==
 
221
      integer, intent(in)           :: g_zdens_ao(nsets)
 
222
C      integer, intent(in)           :: g_movecs_ao_gs(nsets)
 
223
 
 
224
 
 
225
C     == Parameters ==
 
226
      character(*),parameter :: pname="rt_tddft_movecs_zdens: "
 
227
 
 
228
 
 
229
C     == External ==
 
230
      logical, external  :: movecs_read_header, movecs_read
 
231
 
 
232
      
 
233
C     == Variables ==
 
234
      character*256 rt_movecs_fname
 
235
      integer g_movecs_pad          !double, ns_ao x ns_ao w/ 0's for canorg removed lindeps
 
236
      integer dtype, dim1, dim2
 
237
      integer me
 
238
      integer g_densre_ao
 
239
 
 
240
C     (movecs header stuff)
 
241
      character*256 mo_title, mo_basis_name, mo_scftype
 
242
      integer mo_nbf, mo_nsets, mo_ns_mo(nsets)
 
243
 
 
244
      integer locc, iocc, levals, ievals  !MO occupations and eigenvalues (from movecs file)
 
245
 
 
246
      integer g_tmp
 
247
      integer iorb
 
248
      integer is
 
249
      double precision occ
 
250
 
 
251
      me = ga_nodeid ()
 
252
 
 
253
 
 
254
C
 
255
C     Checks
 
256
C
 
257
C      if ((ispin .lt. 1).or.(ispin .gt. nsets))
 
258
C     $     call errquit (pname//"invalid ispin supplied", ispin, 0)
 
259
 
 
260
      if ((nsets .ne. 1).and.(nsets .ne. 2))
 
261
     $     call errquit (pname//"invalid nsets supplied", nsets, 0)
 
262
 
 
263
      if (nsets .ne. 1)
 
264
     $     call errquit (pname//"nsets > 1 needs to be implemented",0,0)
 
265
 
 
266
      do is = 1, nsets
 
267
         call ga_check_handle (g_zdens_ao(is),
 
268
     $        "third argument of "//pname//"not a valid GA")
 
269
         call ga_inquire (g_zdens_ao(is), dtype, dim1, dim2)
 
270
      
 
271
         if (dtype .ne. mt_dcpl) call errquit (pname//
 
272
     $        "expecting complex-valued GA as third argument", 0, 0)
 
273
         if (dim1 .ne. dim2)
 
274
     $        call errquit (pname//"dim1 must equal dim2", 0, 0)
 
275
         if (dim1 .ne. params%ns_ao)
 
276
     $        call errquit (pname//
 
277
     $        "bad size P--expecting ns_ao x ns_ao", 0, 0)
 
278
      enddo
 
279
 
 
280
 
 
281
 
 
282
C
 
283
C     Allocation
 
284
C
 
285
      if (.not. ga_create(mt_dbl,params%ns_ao,params%ns_ao,
 
286
     $     "movecs", 0, 0, g_movecs_pad))
 
287
     $     call errquit ("couldnt create movecs_pad", 0, GA_ERR)
 
288
 
 
289
C
 
290
C     Just an alias to the ouput GA; note padded with 0's if have
 
291
C     lindep)
 
292
C     
 
293
CXXX  [KAL]: ONLY FOR nsets = 1
 
294
C      g_movecs_pad = g_movecs_ao_gs(1)
 
295
 
 
296
 
 
297
      if (.not. ga_duplicate(g_movecs_pad, g_densre_ao, "real P"))
 
298
     $     call errquit ("couldnt duplicate movecs", 0, GA_ERR)
 
299
 
 
300
      if (.not. ga_duplicate(g_movecs_pad, g_tmp, "tmp"))
 
301
     $     call errquit ("couldnt duplicate movecs", 0, GA_ERR)
 
302
 
 
303
 
 
304
C
 
305
C     Read in header to get file name, check, then read in movecs.  Note
 
306
C     that the g_movecs_pad is ns_ao x ns_ao, with the last few columns
 
307
C     possibly 0 (if lindep), which is how the SCF code does it, which
 
308
C     is unlike my way, which has ns_mo x ns_mo.
 
309
C
 
310
C
 
311
      call ga_zero (g_movecs_pad)
 
312
 
 
313
      call rt_tddft_movecs_fname (params, rt_movecs_fname)
 
314
 
 
315
      if (.not. movecs_read_header (rt_movecs_fname, mo_title,
 
316
     $     mo_basis_name, mo_scftype, mo_nbf, mo_nsets,
 
317
     $     mo_ns_mo, nsets))
 
318
     $     call errquit (pname//"Failed to read movecs header", 0, 0)
 
319
 
 
320
 
 
321
C
 
322
C     Check that movecs are legit.
 
323
C
 
324
      if (mo_scftype .ne. "dft")
 
325
     $     call errquit (pname//
 
326
     $     'Initial movecs should have scftype "dft"', 0, 0)
 
327
 
 
328
      if (mo_nbf .ne. params%ns_ao)
 
329
     $     call errquit (pname//
 
330
     $     'Initial movecs wrong size: mo_nbf /= ns_ao', mo_nbf, 0)
 
331
 
 
332
      do is = 1, nsets
 
333
         if (mo_ns_mo(is) .ne. params%ns_mo)
 
334
     $        call errquit (pname//
 
335
     $        'Initial movecs wrong size: mo_ns_mo /= ns_mo',
 
336
     $        mo_ns_mo(is), 0)
 
337
      enddo
 
338
      
 
339
      if (mo_nsets .ne. nsets)
 
340
     $     call errquit (pname//"Wrong number of initial movecs,",
 
341
     $     mo_nsets, 0)
 
342
 
 
343
      
 
344
      call rt_tddft_movecs_print_header (params, rt_movecs_fname,
 
345
     $     mo_title, mo_basis_name, mo_scftype, mo_nbf,
 
346
     $     mo_nsets, mo_ns_mo)
 
347
 
 
348
 
 
349
 
 
350
C
 
351
C     Allocate buffers and movecs (ns_ao x ns_ao padded with zero *not* ns_mo x ns_mo).
 
352
C
 
353
      if (.not.ma_push_get(mt_dbl, params%ns_ao, 'occ', locc, iocc))
 
354
     &     call errquit(pname//'cannot allocate occ',0, MA_ERR)
 
355
 
 
356
      if (.not.ma_push_get(mt_dbl, params%ns_ao, 'evals',
 
357
     $     levals, ievals))
 
358
     &     call errquit(pname//'cannot allocate evals',0, MA_ERR)
 
359
 
 
360
 
 
361
C
 
362
C     Read in movecs (note ispin).
 
363
C
 
364
      if (.not. movecs_read (rt_movecs_fname, 1, dbl_mb(iocc),  !note 1!!!
 
365
     $     dbl_mb(ievals), g_movecs_pad))
 
366
     $     call errquit (pname//"Failed to read movecs data", 0, 0)
 
367
 
 
368
      call rt_tddft_movecs_print_evals (params,
 
369
     $     dbl_mb(iocc), dbl_mb(ievals))
 
370
 
 
371
 
 
372
C
 
373
C     Mask the movecs and multiply into the unmasked transposed movecs
 
374
C     to make the real-valued dens mat in AO basis, then cast to
 
375
C     complex.
 
376
C
 
377
      call ga_zero (g_tmp)
 
378
 
 
379
      do iorb = 1, params%ns_ao
 
380
         occ = dbl_mb(iocc + iorb - 1)
 
381
         call ga_fill_patch (g_tmp, 1, params%ns_ao,
 
382
     $        iorb, iorb, occ)
 
383
      enddo
 
384
 
 
385
      call ga_sync () !XXX needed?
 
386
 
 
387
      call ga_elem_multiply (g_tmp, g_movecs_pad, g_tmp) ! g_tmp now holds masked movecs
 
388
      call ga_zero (g_densre_ao)
 
389
      call ga_dgemm ("N", "T", params%ns_ao, params%ns_ao, params%ns_ao,
 
390
     $     1d0, g_tmp, g_movecs_pad, 0d0, g_densre_ao)
 
391
 
 
392
      call convert_d2z(1d0, g_densre_ao, 0d0, g_densre_ao, g_zdens_ao)
 
393
 
 
394
C      call ga_print (g_movecs_pad)
 
395
 
 
396
 
 
397
C     Clean up
 
398
C
 
399
      if (.not. ga_destroy (g_movecs_pad))
 
400
     $     call errquit (pname//"failed to destroy movecs",0,0)
 
401
 
 
402
      if (.not. ga_destroy (g_densre_ao))
 
403
     $     call errquit (pname//"failed to destroy densao_tmp",0,0)
 
404
 
 
405
      if (.not. ma_pop_stack (levals))
 
406
     $     call errquit (pname//"Failed to pop evals", 0, ma_err)
 
407
      
 
408
      if (.not. ma_pop_stack (locc))
 
409
     $     call errquit (pname//"Failed to pop occ", 0, ma_err)
 
410
 
 
411
      if (.not. ga_destroy (g_tmp))
 
412
     $     call errquit ("couldnt destroy tmp", 0, GA_ERR)
 
413
 
 
414
      end subroutine 
 
415
 
 
416
 
 
417
 
 
418
 
 
419
 
 
420
CXXX  [KAL]: REMOVE ONCE UNIFIED WITH CLOSEDSHELL ROUTINE
 
421
C====================================================================
 
422
C
 
423
C     rt_tddft_os_movecs_zdens.F
 
424
C
 
425
C     Read in initial state movecs from file and convert to OPEN SHELL
 
426
C     complex dens matrix in AO basis.  Although the output dens mat is
 
427
C     complex data type, it is pure real (as the movecs are from the SCF
 
428
C     and thus pure real).
 
429
C
 
430
C
 
431
      subroutine rt_tddft_os_movecs_zdens (params, g_zdens_ao)
 
432
      implicit none
 
433
      
 
434
#include "errquit.fh"
 
435
#include "mafdecls.fh"
 
436
#include "stdio.fh"
 
437
#include "global.fh"
 
438
#include "msgids.fh"
 
439
#include "geom.fh"
 
440
#include "util.fh"
 
441
#include "cdft.fh"
 
442
#include "rtdb.fh"
 
443
#include "rt_tddft.fh"
 
444
#include "matutils.fh"
 
445
 
 
446
      
 
447
C     == Inputs ==
 
448
      type (rt_params_t), intent(in) :: params
 
449
 
 
450
      
 
451
C     == Outputs ==
 
452
      integer, intent(in) :: g_zdens_ao(2)     !dble complex, ns_ao x ns_ao; alpha, beta
 
453
 
 
454
 
 
455
C     == Parameters ==
 
456
      character(*),parameter :: pname="rt_tddft_os_movecs_zdens: "
 
457
 
 
458
      
 
459
C     == External ==
 
460
      logical, external  :: movecs_read_header, movecs_read
 
461
 
 
462
      
 
463
C     == Variables ==
 
464
      character*256 rt_movecs_fname
 
465
      integer dtype, dim1, dim2
 
466
      integer me
 
467
      integer g_movecs_pad
 
468
      integer g_densao_tmp
 
469
      logical ok
 
470
 
 
471
C     (movecs header stuff)
 
472
      character*256 mo_title, mo_basis_name, mo_scftype
 
473
      integer mo_nbf, mo_nsets, mo_ns_mo(2) !2 for alpha, beta spins
 
474
 
 
475
      integer locc, iocc, levals, ievals  !MO occupations and eigenvalues (from movecs file)
 
476
      integer is
 
477
 
 
478
 
 
479
      call rt_tddft_os_confirm (params)
 
480
      
 
481
      me = ga_nodeid ()
 
482
      
 
483
      if (me.eq.0) then
 
484
         write (luout, *) ""
 
485
         write (luout, *) ""
 
486
         call util_print_centered (luout,
 
487
     $        "Initial state: Imported open shell MO eigenvectors",
 
488
     $        20, .true.)
 
489
      endif
 
490
         
 
491
 
 
492
C
 
493
C     Check GAs
 
494
C
 
495
      call ga_check_handle (g_zdens_ao(1),
 
496
     $     "third argument of "//pname//"not a valid GA")
 
497
      call ga_inquire (g_zdens_ao(1), dtype, dim1, dim2)
 
498
      if (dtype .ne. mt_dcpl) call errquit (pname//
 
499
     $     "expecting complex-valued GA as third argument", 0, 0)
 
500
      if (dim1 .ne. dim2)
 
501
     $     call errquit (pname//"dim1 must equal dim2", 0, 0)
 
502
      if (dim1 .ne. params%ns_ao)
 
503
     $     call errquit (pname//"bad size P--expecting ns_ao x ns_ao",
 
504
     $     0, 0)
 
505
 
 
506
      call ga_check_handle (g_zdens_ao(2),
 
507
     $     "third argument of "//pname//"not a valid GA")
 
508
      call ga_inquire (g_zdens_ao(2), dtype, dim1, dim2)
 
509
      if (dtype .ne. mt_dcpl) call errquit (pname//
 
510
     $     "expecting complex-valued GA as third argument", 0, 0)
 
511
      if (dim1 .ne. dim2)
 
512
     $     call errquit (pname//"dim1 must equal dim2", 0, 0)
 
513
      if (dim1 .ne. params%ns_ao)
 
514
     $     call errquit (pname//"bad size P--expecting ns_ao x ns_ao",
 
515
     $     0, 0)
 
516
 
 
517
 
 
518
C
 
519
C     Allocation
 
520
C
 
521
      if (.not. ga_create(mt_dbl,params%ns_ao,params%ns_ao,
 
522
     $     "movecs", 0, 0, g_movecs_pad))
 
523
     $     call errquit ("couldnt create movecs_pad", 0, GA_ERR)
 
524
 
 
525
      if (.not. ga_duplicate(g_movecs_pad, g_densao_tmp, "real P"))
 
526
     $     call errquit ("couldnt duplicate movecs", 0, GA_ERR)
 
527
 
 
528
      
 
529
 
 
530
C
 
531
C     Read in header to get file name, check, then read in movecs.  Note
 
532
C     that the g_movecs_pad_pad is ns_ao x ns_ao, with the last few columns
 
533
C     possibly 0 (if lindep), which is how the SCF code does it, which
 
534
C     is unlike my way, which has ns_mo x ns_mo.  The 2 is for open shell.
 
535
C
 
536
      call rt_tddft_movecs_fname (params, rt_movecs_fname)
 
537
 
 
538
 
 
539
      if (.not. movecs_read_header (rt_movecs_fname, mo_title,
 
540
     $     mo_basis_name, mo_scftype, mo_nbf, mo_nsets,
 
541
     $     mo_ns_mo, 2))
 
542
     $     call errquit (pname//"Failed to read movecs header", 0, 0)
 
543
 
 
544
C      call rt_tddft_movecs_print_header (params, rt_movecs_fname,
 
545
C     $     mo_title, mo_basis_name, mo_scftype, mo_nbf)
 
546
 
 
547
 
 
548
C
 
549
C     Check that movecs are legit.
 
550
C
 
551
      if (mo_scftype .ne. "dft")
 
552
     $     call errquit (pname//
 
553
     $     'Initial movecs should have scftype "dft"', 0, 0)
 
554
 
 
555
      if (mo_nbf .ne. params%ns_ao)
 
556
     $     call errquit (pname//
 
557
     $     'Initial movecs wrong size: mo_nbf /= ns_ao', 0, 0)
 
558
 
 
559
      if (mo_ns_mo(1) .ne. params%ns_mo)
 
560
     $     call errquit (pname//
 
561
     $     'Initial movecs wrong size: mo_ns_mo(1) /= ns_mo', 0, 0)
 
562
 
 
563
      if (mo_ns_mo(2) .ne. params%ns_mo)
 
564
     $     call errquit (pname//
 
565
     $     'Initial movecs wrong size: mo_ns_mo(2) /= ns_mo', 0, 0)
 
566
      
 
567
C      call rt_tddft_print_warning ("Didnt check mo_ns_mo(2) == ns_mo")
 
568
 
 
569
      if (mo_nsets .ne. 2)
 
570
     $     call errquit (pname//"Wrong number of initial movecs,"//
 
571
     $     " mo_nsets should be 2 for open shell", 0, 0)
 
572
 
 
573
 
 
574
C
 
575
C     Allocate tmp buffers for occupations and evals (ns_ao x ns_ao
 
576
C     padded with zero *not* ns_mo x ns_mo).
 
577
C
 
578
      if (.not.ma_push_get(mt_dbl, params%ns_ao, 'occ', locc, iocc))
 
579
     &     call errquit(pname//'cannot allocate occ 1',0, MA_ERR)
 
580
 
 
581
      if (.not.ma_push_get(mt_dbl, params%ns_ao, 'evals',
 
582
     $     levals, ievals))
 
583
     &     call errquit(pname//'cannot allocate evals',0, MA_ERR)
 
584
 
 
585
 
 
586
C
 
587
C     Read in movecs.  Note we loop over the two spins since closed
 
588
C     shell, and we recycle all the MA and GA space, and finally store
 
589
C     in g_zdens_ao which has two components (alpha, beta spins).
 
590
C
 
591
      call ga_zero (g_zdens_ao(1))
 
592
      call ga_zero (g_zdens_ao(2))
 
593
 
 
594
C     (alpha part)
 
595
      call ga_zero (g_movecs_pad)
 
596
         
 
597
      if (.not. movecs_read (rt_movecs_fname, 1, dbl_mb(iocc),
 
598
     $     dbl_mb(ievals), g_movecs_pad))
 
599
     $     call errquit (pname//"Failed to read movecs data", 0, 0)
 
600
      
 
601
      call rt_tddft_movecs_print_evals (params,
 
602
     $     dbl_mb(iocc), dbl_mb(ievals))
 
603
      
 
604
      call ga_zero (g_densao_tmp)
 
605
      call ga_dgemm ("N", "T", params%ns_ao, params%ns_ao,
 
606
     $     params%nalpha, 1d0, g_movecs_pad, g_movecs_pad,
 
607
     $     0d0, g_densao_tmp)
 
608
      
 
609
      call convert_d2z(1d0, g_densao_tmp, 0d0, g_densao_tmp,
 
610
     $     g_zdens_ao(1))
 
611
 
 
612
 
 
613
C     (beta part)
 
614
      call ga_zero (g_movecs_pad)
 
615
         
 
616
      if (.not. movecs_read (rt_movecs_fname, 2, dbl_mb(iocc),
 
617
     $     dbl_mb(ievals), g_movecs_pad))
 
618
     $     call errquit (pname//"Failed to read movecs data", 0, 0)
 
619
      
 
620
      call rt_tddft_movecs_print_evals (params,
 
621
     $     dbl_mb(iocc), dbl_mb(ievals))
 
622
      
 
623
      
 
624
      call ga_zero (g_densao_tmp)
 
625
      call ga_dgemm ("N", "T", params%ns_ao, params%ns_ao,
 
626
     $     params%nbeta, 1d0, g_movecs_pad, g_movecs_pad,
 
627
     $     0d0, g_densao_tmp)
 
628
      
 
629
      call convert_d2z(1d0, g_densao_tmp, 0d0, g_densao_tmp,
 
630
     $     g_zdens_ao(2))
 
631
      
 
632
 
 
633
C
 
634
C     Clean up
 
635
C
 
636
      if (.not. ga_destroy (g_movecs_pad))
 
637
     $     call errquit (pname//"failed to destroy movecs",0,0)
 
638
 
 
639
      if (.not. ga_destroy (g_densao_tmp))
 
640
     $     call errquit (pname//"failed to destroy densao_tmp",0,0)
 
641
 
 
642
      if (.not. ma_pop_stack (levals))
 
643
     $     call errquit (pname//"Failed to pop evals", 0, ma_err)
 
644
      
 
645
      if (.not. ma_pop_stack (locc))
 
646
     $     call errquit (pname//"Failed to pop occ", 0, ma_err)
 
647
 
 
648
 
 
649
      end subroutine