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

« back to all changes in this revision

Viewing changes to src/nwdft/zora/dft_zora_EPR-NMR_tools.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
* 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:
4
4
     &                           nbf,      ! in  : nr basis functions
5
5
     &                           geom,     ! in  : geometry handle
6
6
     &                           mcenters) ! in  : nr. atoms
 
7
c
 
8
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
9
 
7
10
      implicit none
8
11
 
9
12
#include "rtdb.fh" 
60
63
      subroutine get_chi_centers_ga(g_chi_cntr, ! out
61
64
     &                              basis,      ! in  : basis handle
62
65
     &                              nbf,geom,mcenters)
 
66
c
 
67
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
68
 
63
69
      implicit none
64
70
 
65
71
#include "rtdb.fh" 
142
148
     &                         ind_b,   ! from kab=123,231,312
143
149
     &                         g_tmp2,  ! scratch
144
150
     &                         g_N_scld)! output
 
151
c
 
152
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
153
 
145
154
       implicit none
146
155
#include "errquit.fh"
147
156
#include "mafdecls.fh"
164
173
      subroutine get_scld_A(g_A,  ! ga-arr to scale - OUT
165
174
     &                      g_R,  ! scaling arr
166
175
     &                      g_tmp)! scratch arr 
 
176
c
 
177
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
178
 
167
179
       implicit none
168
180
#include "errquit.fh"
169
181
#include "mafdecls.fh"
199
211
     &            g_para1, ! in: para1 tensor
200
212
     &              g_h01, ! in: h01 AO matrix
201
213
     &              g_Fji) ! in: Perturbed Fock matrix
 
214
c
 
215
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
216
 
202
217
      implicit none
203
218
#include "errquit.fh"
204
219
#include "mafdecls.fh"
371
386
     &            g_para1, ! out: para1 tensor
372
387
     &              g_h01, ! out: h01 AO matrix
373
388
     &              g_Fji) ! out: Perturbed Fock matrix
 
389
c
 
390
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
391
 
374
392
      implicit none
375
393
#include "errquit.fh"
376
394
#include "global.fh"
597
615
     &            g_AtNr1, ! in: list of atoms to calc. shieldings
598
616
     &              g_dia, ! in: dia A,B tensor
599
617
     &            g_para1) ! in: par A,B tensor
 
618
c
 
619
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
620
 
600
621
      implicit none
601
622
#include "errquit.fh"
602
623
#include "mafdecls.fh"
725
746
     &            g_AtNr1, ! out: list of atoms to calc. shieldings
726
747
     &              g_dia, ! out: dia-A,B   tensor
727
748
     &            g_para1) ! out: par-A,B   tensor    
 
749
c
 
750
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
751
 
728
752
      implicit none
729
753
#include "errquit.fh"
730
754
#include "global.fh"
893
917
     &           vectors,  ! in: MOs
894
918
     &           g_rhs0,   ! in: (ntot,3)       GA matrix
895
919
     &           g_rhs)    ! in: (nocc*nvirt,3) GA matrix
 
920
c
 
921
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
922
 
896
923
      implicit none
897
924
#include "errquit.fh"
898
925
#include "mafdecls.fh"
1009
1036
      dft_zoraCPHF_write = (ok .eq. 1)
1010
1037
      if (ga_nodeid() .eq. 0) then
1011
1038
         write(6,22) filename(1:inp_strlen(filename))
1012
 
 22      format(/' Wrote ZORA CPHF data to ',a/)
 
1039
 22      format(/' Wrote CPHF data to ',a/)
1013
1040
         call util_flush(luout)
1014
1041
      endif
1015
1042
      return
1040
1067
     &           vectors,  ! out: MOs
1041
1068
     &           g_rhs0,   ! out: (ntot,3)       GA matrix
1042
1069
     &           g_rhs)    ! out: (nocc*nvirt,3) GA matrix
 
1070
c
 
1071
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
1072
 
1043
1073
      implicit none
1044
1074
#include "errquit.fh"
1045
1075
#include "global.fh"
1186
1216
      goto 10
1187
1217
      end
1188
1218
c =========== READ/WRITE CPHF (g_rhs) data ==========END
1189
 
c $Id: dft_zora_EPR-NMR_tools.F 21176 2011-10-10 06:35:49Z d3y133 $
 
1219
c $Id: dft_zora_EPR-NMR_tools.F 23379 2013-01-05 23:41:27Z niri $
 
1220
c =========== READ/WRITE CPHF-1 (g_rhs) data ==========START
 
1221
c To be used in aoresponse module: fiao_f1_movecs
 
1222
      logical function dft_CPHF1_write(
 
1223
     &           filename, ! in: filename
 
1224
     &           npol,     ! in: nr polarization
 
1225
     &           nocc,     ! in: nr occupied MOs
 
1226
     &           nvirt,    ! in: nr virtual  MOs
 
1227
     &           ncomp,    ! in: nr. components
 
1228
     &           g_rhs_re, ! in: (nocc*nvirt,3)(ipm) GA matrix
 
1229
     &           g_rhs_im, ! in: (nocc*nvirt,3)(ipm) GA matrix
 
1230
     &           lifetime) ! in: =T if (RE,IM) =F if RE
 
1231
c
 
1232
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
1233
c
 
1234
c Note.- nmo ne nbf if linear dependencies appear.
 
1235
 
 
1236
      implicit none
 
1237
#include "errquit.fh"
 
1238
#include "mafdecls.fh"
 
1239
#include "global.fh"
 
1240
#include "tcgmsg.fh"
 
1241
#include "msgtypesf.h"
 
1242
#include "inp.fh"
 
1243
#include "msgids.fh"
 
1244
#include "cscfps.fh"
 
1245
#include "util.fh"
 
1246
#include "stdio.fh"
 
1247
      character*(*) filename    ! [input] File to write to
 
1248
      integer npol,
 
1249
     &        nocc(npol),nvirt(npol),
 
1250
     &        ispin,ntot,
 
1251
     &        ipm,ncomp,
 
1252
     &        g_rhs_re(ncomp),
 
1253
     &        g_rhs_im(ncomp)
 
1254
      integer unitno
 
1255
      parameter (unitno = 77)
 
1256
      integer l_rhs0,k_rhs0,
 
1257
     &        l_rhs,k_rhs
 
1258
      integer ok,i,j
 
1259
      integer inntsize
 
1260
      integer nrhs,nxyz
 
1261
      logical lifetime
 
1262
 
 
1263
      nxyz=3 ! =x,y,z
 
1264
      inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
 
1265
      call ga_sync()
 
1266
      ok = 0
 
1267
c     Read routines should be consistent with this
 
1268
c     Write out the atomic zora corrections
 
1269
      if (ga_nodeid() .eq. 0) then
 
1270
c     Open the file
 
1271
       open(unitno, status='unknown', form='unformatted',
 
1272
     $      file=filename, err=1000)
 
1273
c     Write out the number of sets and basis functions
 
1274
       write(unitno, err=1001) npol
 
1275
       do i=1,npol
 
1276
        write(unitno, err=1001) nocc(i)
 
1277
       enddo
 
1278
       do i=1,npol
 
1279
        write(unitno, err=1001) nvirt(i)
 
1280
       enddo
 
1281
       write(unitno, err=1001) nxyz
 
1282
c     Allocate the temporary buffer
 
1283
       ntot=0
 
1284
       do ispin=1,npol
 
1285
        ntot=ntot+nocc(ispin)*nvirt(ispin)
 
1286
       enddo
 
1287
       write(unitno, err=1001) ntot
 
1288
       if (.not. ma_alloc_get(mt_dbl,ntot,
 
1289
     &           'dft_CPHF_write',l_rhs,k_rhs))
 
1290
     $  call errquit('dft_CPHF_write: ma failed', 
 
1291
     &               ntot, MA_ERR)
 
1292
       do ipm=1,ncomp
 
1293
        do i=1,2*nxyz  ! write (g_b,g_rhs_sol)
 
1294
         call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs),1)
 
1295
         call ga_get(g_rhs_re(ipm),1,ntot,i,i,dbl_mb(k_rhs),1)
 
1296
         call swrite(unitno,dbl_mb(k_rhs),ntot)
 
1297
        enddo ! end-loop-i
 
1298
        if (lifetime) then
 
1299
        do i=1,2*nxyz ! write (g_b,g_rhs_sol)
 
1300
         call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs),1)
 
1301
         call ga_get(g_rhs_im(ipm),1,ntot,i,i,dbl_mb(k_rhs),1)
 
1302
         call swrite(unitno,dbl_mb(k_rhs),ntot)
 
1303
        enddo ! end-loop-i
 
1304
        endif ! end-if-lifetime
 
1305
       enddo ! end-loop-ipm
 
1306
c
 
1307
c     Deallocate the temporary buffer
 
1308
c ----- Using ma_free_heap ------------ START
 
1309
      if (.not. ma_free_heap(l_rhs))
 
1310
     $  call errquit('dft_CPHF_write: ma free_heap failed', 
 
1311
     &               911, MA_ERR)
 
1312
c ----- Using ma_free_heap ------------ END
 
1313
c     Close the file
 
1314
      close(unitno,err=1002)
 
1315
      ok = 1
 
1316
      end if
 
1317
c     Broadcast status to other nodes
 
1318
 10   call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status
 
1319
      call ga_sync()
 
1320
      dft_CPHF1_write = (ok .eq. 1)
 
1321
      if (ga_nodeid() .eq. 0) then
 
1322
         write(6,22) filename(1:inp_strlen(filename))
 
1323
 22      format(/' dft_CPHF_write: Wrote aoresponse g_rhs data to ',a/)
 
1324
         call util_flush(luout)
 
1325
      endif
 
1326
      return
 
1327
 1000 write(6,*) 'dft_CPHF_write: failed to open ',
 
1328
     $     filename(1:inp_strlen(filename))
 
1329
      call util_flush(luout)
 
1330
      ok = 0
 
1331
      goto 10
 
1332
 1001 write(6,*) 'dft_CPHF_write: failed to write ',
 
1333
     $     filename(1:inp_strlen(filename))
 
1334
      call util_flush(luout)
 
1335
      ok = 0
 
1336
      close(unitno,err=1002)
 
1337
      goto 10
 
1338
 1002 write(6,*) 'dft_CPHF_write: failed to close',
 
1339
     $     filename(1:inp_strlen(filename))
 
1340
      call util_flush(luout)
 
1341
      ok = 0
 
1342
      goto 10
 
1343
      end
 
1344
 
 
1345
      logical function dft_CPHF1_read(
 
1346
     &           filename, !  in: filename
 
1347
     &           npol,     !  in: nr polarization
 
1348
     &           nocc,     !  in: nr occupied MOs
 
1349
     &           nvirt,    !  in: nr virtual  MOs
 
1350
     &           ncomp,    ! out: nr. components
 
1351
     &           g_rhs_re, ! out: (nocc*nvirt,3)(ipm) GA matrix
 
1352
     &           g_rhs_im, ! out: (nocc*nvirt,3)(ipm) GA matrix
 
1353
     &           lifetime) ! out: =T if (RE,IM) =F if RE
 
1354
c
 
1355
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
1356
 
 
1357
      implicit none
 
1358
#include "errquit.fh"
 
1359
#include "global.fh"
 
1360
#include "tcgmsg.fh"
 
1361
#include "msgtypesf.h"
 
1362
#include "mafdecls.fh"
 
1363
#include "msgids.fh"
 
1364
#include "cscfps.fh"
 
1365
#include "inp.fh"
 
1366
#include "util.fh"
 
1367
#include "stdio.fh"
 
1368
      character*(*) filename    ! [input] File to write to
 
1369
      integer npol,
 
1370
     &        nocc(npol),nvirt(npol),
 
1371
     &        ispin,ntot,
 
1372
     &        ipm,ncomp,
 
1373
     &        g_rhs_re(ncomp),
 
1374
     &        g_rhs_im(ncomp)
 
1375
      integer unitno
 
1376
      parameter (unitno = 77)
 
1377
      integer l_rhs,k_rhs
 
1378
      integer ok,i,j
 
1379
      integer inntsize
 
1380
      integer nrhs,nxyz
 
1381
      integer npol_read,nxyz_read,ntot_read,
 
1382
     &        nocc_read(2),nvirt_read(2)
 
1383
      logical lifetime
 
1384
      nxyz=3 ! =x,y,z   
 
1385
c     Initialise to invalid MA handle
 
1386
      inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
 
1387
      call ga_sync()
 
1388
      ok = 0
 
1389
      if (ga_nodeid() .eq. 0) then
 
1390
c      Print a message indicating the file being read
 
1391
       write(6,22) filename(1:inp_strlen(filename))
 
1392
 22    format(/' dft_CPHF_read:Read aoresponse g_rhs data from ',a/)
 
1393
       call util_flush(luout)
 
1394
c      Open the file
 
1395
       open(unitno, status='old', form='unformatted', file=filename,
 
1396
     $        err=1000)
 
1397
c      Read in some basics to check if they are consistent with the calculation
 
1398
       read(unitno, err=1001, end=1001) npol_read
 
1399
       do i=1,npol_read
 
1400
        read(unitno, err=1001, end=1001) nocc_read(i)
 
1401
       enddo
 
1402
       do i=1,npol_read
 
1403
        read(unitno, err=1001, end=1001) nvirt_read(i)
 
1404
       enddo
 
1405
       read(unitno, err=1001, end=1001) nxyz_read
 
1406
c      Error checks
 
1407
       if ((nxyz_read  .ne. nxyz) .or.
 
1408
     &     (npol_read  .ne. npol)) goto 1003
 
1409
       ntot=0
 
1410
       do ispin=1,npol
 
1411
        ntot=ntot+nocc(ispin)*nvirt(ispin)
 
1412
       enddo
 
1413
       read(unitno, err=1001, end=1001) ntot_read
 
1414
       if (.not. ma_alloc_get(mt_dbl,ntot,
 
1415
     &           'dft_CPHF_read',l_rhs,k_rhs))
 
1416
     $  call errquit('dft_CPHF_read: ma failed', 
 
1417
     &               ntot, MA_ERR)
 
1418
 
 
1419
       do ipm=1,ncomp
 
1420
        do i=1,nxyz ! skip 1st subspace
 
1421
         call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs),1)
 
1422
         call sread(unitno,dbl_mb(k_rhs),ntot)
 
1423
        enddo 
 
1424
        do i=nxyz+1,2*nxyz ! read 2nd subspace and copy to g_rhs_re
 
1425
         call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs),1)
 
1426
         call sread(unitno,dbl_mb(k_rhs),ntot)
 
1427
         call ga_put(g_rhs_re(ipm),1,ntot,i,i,dbl_mb(k_rhs),1)
 
1428
        enddo ! end-loop-i
 
1429
        if (lifetime) then
 
1430
        do i=1,nxyz ! skip 1st subspace
 
1431
         call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs),1)
 
1432
         call sread(unitno,dbl_mb(k_rhs),ntot)
 
1433
        enddo ! end-loop-i
 
1434
        do i=nxyz+1,2*nxyz ! read 2nd subspace and copy to g_rhs_im
 
1435
         call dcopy(ntot,0.0d0,0,dbl_mb(k_rhs),1)
 
1436
         call sread(unitno,dbl_mb(k_rhs),ntot)
 
1437
         call ga_put(g_rhs_im(ipm),1,ntot,i,i,dbl_mb(k_rhs),1)
 
1438
        enddo ! end-loop-i
 
1439
        endif ! end-if-lifetime
 
1440
       enddo ! end-loop-ipm
 
1441
c
 
1442
c     Deallocate the temporary buffer
 
1443
c ----- Using ma_free_heap ------------ START
 
1444
       if (.not. ma_free_heap(l_rhs))
 
1445
     $  call errquit('dft_CPHF_read: ma free_heap failed', 
 
1446
     &               911, MA_ERR)
 
1447
c ----- Using ma_free_heap ------------ END
 
1448
c      Close the file
 
1449
       close(unitno,err=1002)
 
1450
       ok = 1
 
1451
      end if
 
1452
c
 
1453
c     Broadcast status to other nodes
 
1454
 10   call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status
 
1455
      call ga_sync()
 
1456
      dft_CPHF1_read = ok .eq. 1
 
1457
      return
 
1458
 1000 write(6,*) 'dft_CPHF_read: failed to open',
 
1459
     $     filename(1:inp_strlen(filename))
 
1460
      call util_flush(luout)
 
1461
      ok = 0
 
1462
      goto 10
 
1463
 1001 write(6,*) 'dft_CPHF_read: failed to read',
 
1464
     $     filename(1:inp_strlen(filename))
 
1465
      call util_flush(luout)
 
1466
      ok = 0
 
1467
      close(unitno,err=1002)
 
1468
      goto 10
 
1469
 1003 write(6,*) 
 
1470
     & 'dft_CPHF_read: file inconsistent with calculation',
 
1471
     $     filename(1:inp_strlen(filename))
 
1472
      call util_flush(luout)
 
1473
      ok = 0
 
1474
      close(unitno,err=1002)
 
1475
      goto 10
 
1476
 1002 write(6,*) 'dft_CPHF_read: failed to close',
 
1477
     $     filename(1:inp_strlen(filename))
 
1478
      call util_flush(luout)
 
1479
      ok = 0
 
1480
      goto 10
 
1481
      end
 
1482
c 000000000000000000000000000000000000000000000000000000
 
1483
c 000000000000000000000000000000000000000000000000000000
 
1484
      logical function dft_CPHF2_write(
 
1485
     &           filename, ! in: filename
 
1486
     &           n,        ! in: sum_{i=1,npol} nocc(i)*nvirt(i)
 
1487
     &           ncomp,    ! in: nr. components
 
1488
     &           nvec,     ! in: nr. of directions = 3
 
1489
     &           n1,       ! in: =n*ncomp
 
1490
     &           nsub,     ! in: last subspace index (nsub+1)= nr of subspaces stored
 
1491
     &           nsub_file,! ub: subspace counter
 
1492
     &           g_z1,     ! in: history matrix z
 
1493
     &           g_Az1)    ! in: history matrix Az
 
1494
c
 
1495
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
1496
c
 
1497
c Note.- nsub = cc * nvec (nvec=3=x,y,z) multiple of nvec
 
1498
c        dim(g_z1) =(n1,maxsub)   maxsub=nvec*11 default
 
1499
c        dim(g_Az1)=(n1,maxsub)   maxsub=nvec*11 default
 
1500
c        nsub <= maxsub
 
1501
c -->It will write only subscape nsub for (z1,Az1)
 
1502
      implicit none
 
1503
#include "errquit.fh"
 
1504
#include "mafdecls.fh"
 
1505
#include "global.fh"
 
1506
#include "tcgmsg.fh"
 
1507
#include "msgtypesf.h"
 
1508
#include "inp.fh"
 
1509
#include "msgids.fh"
 
1510
#include "cscfps.fh"
 
1511
#include "util.fh"
 
1512
#include "stdio.fh"
 
1513
      character*(*) filename    ! [input] File to write to
 
1514
      character*255 filename_mini ! only to store nblocks
 
1515
      integer g_z1,g_Az1
 
1516
      integer unitno
 
1517
      parameter (unitno = 77)
 
1518
      integer l_dat,k_dat,
 
1519
     &        l_z,k_z
 
1520
      integer ok,i,j,nblock,nblock_file,idat
 
1521
      integer inntsize,g_xre,g_xim,m1,iter
 
1522
      integer n,ncomp,nvec,n1,nsub,nsub_file
 
1523
      double precision val_re,val_im
 
1524
      external conv2reim_z1 ! defined in ga_lkain_2cpl3.F
 
1525
 
 
1526
      inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
 
1527
      call ga_sync()
 
1528
      ok = 0
 
1529
c     Read routines should be consistent with this
 
1530
c     Write out the atomic zora corrections
 
1531
      if (ga_nodeid() .eq. 0) then
 
1532
c +++++++++ store nsub +++++++++ START
 
1533
      write(filename_mini,30) trim(filename),'_nblock'
 
1534
 30   format(a,4a)
 
1535
c      write(*,*) 'Writing mini:',filename_mini
 
1536
      nblock=nsub/3+1
 
1537
      nblock_file=nsub_file/3+1
 
1538
c      write(*,2) nblock,nblock_file
 
1539
c 2    format('Writing (nblock,nblock_file)=(',i5,',',i5,')')
 
1540
       open(unitno, status='unknown', form='unformatted',
 
1541
     $      file=filename_mini, err=1000)
 
1542
       write(unitno, err=1001) nblock_file
 
1543
       close(unitno,err=1002)  
 
1544
c +++++++++ store nsub +++++++++ END
 
1545
c     Open the file
 
1546
       open(unitno, status='unknown', form='unformatted',
 
1547
     $      file=filename, err=1000,position='append')
 
1548
c     Write out the number of sets and basis functions
 
1549
       write(unitno, err=1001) n
 
1550
       write(unitno, err=1001) ncomp
 
1551
       write(unitno, err=1001) nvec
 
1552
       write(unitno, err=1001) n1
 
1553
       write(unitno, err=1001) nsub_file
 
1554
c     Allocate the temporary buffer
 
1555
       if (.not. ma_alloc_get(mt_dcpl,n1,
 
1556
     &           'dft_CPHF2_write',l_z,k_z))
 
1557
     $  call errquit('dft_CPHF2_write: ma failed', 
 
1558
     &               n1, MA_ERR)
 
1559
       if (.not. ma_alloc_get(mt_dbl,n1,
 
1560
     &           'dft_CPHF2_write',l_dat,k_dat))
 
1561
     $  call errquit('dft_CPHF2_write: ma failed', 
 
1562
     &               n1, MA_ERR)
 
1563
c ------- write g_z1 ---------------------- START
 
1564
        do i=1,nvec 
 
1565
         m1=(nblock-1)*nvec+i
 
1566
         call ga_get(g_z1,1,n1,m1,m1,dcpl_mb(k_z),1)
 
1567
         do idat=1,n1
 
1568
          val_re=dreal(dcpl_mb(k_z+idat-1))
 
1569
          dbl_mb(k_dat+idat-1)=val_re
 
1570
         enddo ! end-loop-idat
 
1571
         call swrite(unitno,dbl_mb(k_dat),n1)
 
1572
         do idat=1,n1
 
1573
          val_im=dimag(dcpl_mb(k_z+idat-1))
 
1574
          dbl_mb(k_dat+idat-1)=val_im
 
1575
         enddo ! end-loop-idat
 
1576
         call swrite(unitno,dbl_mb(k_dat),n1)
 
1577
        enddo ! end-loop-i
 
1578
c ------- write g_z1 ---------------------- END 
 
1579
c ------- write g_Az1 ---------------------- START
 
1580
        do i=1,nvec 
 
1581
         m1=(nblock-1)*nvec+i
 
1582
         call ga_get(g_Az1,1,n1,m1,m1,dcpl_mb(k_z),1)
 
1583
         do idat=1,n1
 
1584
          val_re=dreal(dcpl_mb(k_z+idat-1))
 
1585
          dbl_mb(k_dat+idat-1)=val_re
 
1586
         enddo ! end-loop-idat
 
1587
         call swrite(unitno,dbl_mb(k_dat),n1)
 
1588
         do idat=1,n1
 
1589
          val_im=dimag(dcpl_mb(k_z+idat-1))
 
1590
          dbl_mb(k_dat+idat-1)=val_im
 
1591
         enddo ! end-loop-idat
 
1592
         call swrite(unitno,dbl_mb(k_dat),n1)
 
1593
        enddo ! end-loop-i
 
1594
c ------- write g_Az1 ---------------------- END
 
1595
c     Deallocate the temporary buffer
 
1596
c ----- Using ma_free_heap ------------ START
 
1597
      if (.not. ma_free_heap(l_dat))
 
1598
     $  call errquit('dft_CPHF2_write: ma free_heap failed', 
 
1599
     &               911, MA_ERR)
 
1600
      if (.not. ma_free_heap(l_z))
 
1601
     $  call errquit('dft_CPHF2_write: ma free_heap failed', 
 
1602
     &               911, MA_ERR)
 
1603
c ----- Using ma_free_heap ------------ END
 
1604
c     Close the file
 
1605
      close(unitno,err=1002)
 
1606
      ok = 1
 
1607
      end if
 
1608
c     Broadcast status to other nodes
 
1609
 10   call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status
 
1610
      call ga_sync()
 
1611
      dft_CPHF2_write = (ok .eq. 1)
 
1612
      if (ga_nodeid() .eq. 0) then
 
1613
        iter=nsub/3+1 ! estimate iter from nsub
 
1614
         write(6,22) filename(1:inp_strlen(filename)),
 
1615
     &               iter,nsub,nsub_file             
 
1616
 22      format(' dft_CPHF2_write: Wrote ',a,
 
1617
     &           ' (iter,nsub,nsub_file)=('
 
1618
     &          i4,',',i5,',',i5,')')
 
1619
         call util_flush(luout)
 
1620
      endif
 
1621
      return
 
1622
 1000 write(6,*) 'dft_CPHF2_write: failed to open ',
 
1623
     $     filename(1:inp_strlen(filename))
 
1624
      call util_flush(luout)
 
1625
      ok = 0
 
1626
      goto 10
 
1627
 1001 write(6,*) 'dft_CPHF2_write: failed to write ',
 
1628
     $     filename(1:inp_strlen(filename))
 
1629
      call util_flush(luout)
 
1630
      ok = 0
 
1631
      close(unitno,err=1002)
 
1632
      goto 10
 
1633
 1002 write(6,*) 'dft_CPHF2_write: failed to close',
 
1634
     $     filename(1:inp_strlen(filename))
 
1635
      call util_flush(luout)
 
1636
      ok = 0
 
1637
      goto 10
 
1638
      end
 
1639
 
 
1640
      logical function dft_CPHF2_read(
 
1641
     &           filename, ! in: filename
 
1642
     &           n,        ! in: sum_{i=1,npol} nocc(i)*nvirt(i)
 
1643
     &           ncomp,    ! in: nr. components
 
1644
     &           nvec,     ! in: nr. of directions = 3
 
1645
     &           n1,       ! in: =n*ncomp
 
1646
     &           nsub,     ! ou: last subspace index (nsub+1)= nr of subspaces stored
 
1647
     &           nsub_read,! ou: last subspace read from file
 
1648
     &           maxsub,   ! in: max subspace of (g_z1,g_Az1)
 
1649
     &           g_z1,     ! ou: history matrix z
 
1650
     &           g_Az1)    ! ou: history matrix Az
 
1651
c
 
1652
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
1653
 
 
1654
      implicit none
 
1655
#include "errquit.fh"
 
1656
#include "global.fh"
 
1657
#include "tcgmsg.fh"
 
1658
#include "msgtypesf.h"
 
1659
#include "mafdecls.fh"
 
1660
#include "msgids.fh"
 
1661
#include "cscfps.fh"
 
1662
#include "inp.fh"
 
1663
#include "util.fh"
 
1664
#include "stdio.fh"
 
1665
      character*(*) filename    ! [input] File to write to
 
1666
      character*255 filename_mini ! only to store nsub
 
1667
      integer ivec,idat,
 
1668
     &        g_z1,g_Az1
 
1669
      integer unitno
 
1670
      parameter (unitno = 77)
 
1671
      integer l_zre,k_zre,
 
1672
     &        l_zim,k_zim
 
1673
      integer ok,i,j,m1,m2,
 
1674
     &        imin,imax,iskip,nskip
 
1675
      integer inntsize
 
1676
      integer n,ncomp,nvec,n1,nsub,maxsub,iblock,
 
1677
     &        n_read,n1_read,nvec_red,ncomp_read,
 
1678
     &        nsub_read,nvec_read,
 
1679
     &        nblocks_ga,nblocks_file
 
1680
      double precision val_re,val_im
 
1681
      double complex val_cmplx
 
1682
      external conv2reim_z1 ! defined in ga_lkain_2cpl3.F
 
1683
 
 
1684
      inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
 
1685
      call ga_sync()
 
1686
      ok = 0
 
1687
      imin=0
 
1688
      imax=0
 
1689
      if (ga_nodeid() .eq. 0) then ! ----- if-ga_nodeid-eq-0 ----- START
 
1690
c +++++++++ read nsub +++++++++ START
 
1691
      write(filename_mini,30) trim(filename),'_nblock'
 
1692
 30   format(a,4a)
 
1693
c      write(*,*) 'Reading mini:',filename_mini
 
1694
       open(unitno, status='old', form='unformatted',
 
1695
     $      file=filename_mini, err=1000)
 
1696
       read(unitno, err=1001, end=1001) nblocks_file
 
1697
       close(unitno,err=1002)  
 
1698
       nblocks_ga=maxsub/nvec
 
1699
       write(*,2) nblocks_file,nblocks_ga
 
1700
 2     format('(nblocks_file,nblocks_ga)=(',i4,',',i4,')')
 
1701
c --- Compare (nblocks_ga,nblocks_file)
 
1702
       if (nblocks_file .le. nblocks_ga-1) then
 
1703
        imin=1
 
1704
        imax=nblocks_file
 
1705
        nskip=0
 
1706
       else
 
1707
        imin=1
 
1708
        imax=nblocks_ga-1
 
1709
        nskip=nblocks_file-(nblocks_ga-1)
 
1710
       endif
 
1711
       write(*,3) imin,imax,nskip
 
1712
 3     format('(imin,imax,nskip)=(',i4,',',i4,',',i4,')')
 
1713
c +++++++++ read nsub +++++++++ END
 
1714
c      Print a message indicating the file being read
 
1715
       write(6,22) filename(1:inp_strlen(filename))
 
1716
 22    format(/' dft_CPHF_read:Read aoresponse g_rhs data from ',a/)
 
1717
       call util_flush(luout)
 
1718
c      Open the file
 
1719
       open(unitno, status='old', form='unformatted', file=filename,
 
1720
     $        err=1000)
 
1721
c      Read in some basics to check if they are consistent with the calculation
 
1722
      if (.not.MA_Push_Get(mt_dbl,n1,'hessv jfacs',l_zre,k_zre))
 
1723
     &     call errquit('conv2complex: cannot allocate zre',
 
1724
     &                  n1, MA_ERR)
 
1725
      if (.not.MA_Push_Get(mt_dbl,n1,'hessv kfacs',l_zim,k_zim))
 
1726
     &     call errquit('conv2complex: cannot allocate zim',
 
1727
     &                  n1, MA_ERR)
 
1728
c ======= skip blocks ======================= START
 
1729
       do iskip=1,nskip
 
1730
        read(unitno, err=1001, end=1001) n_read
 
1731
        read(unitno, err=1001, end=1001) ncomp_read
 
1732
        read(unitno, err=1001, end=1001) nvec_read
 
1733
        read(unitno, err=1001, end=1001) n1_read
 
1734
        read(unitno, err=1001, end=1001) nsub_read
 
1735
c ------- skip g_z1 ---------------------- START
 
1736
         do ivec=1,nvec 
 
1737
          call sread(unitno,dbl_mb(k_zre),n1)
 
1738
          call sread(unitno,dbl_mb(k_zim),n1)
 
1739
         enddo ! end-loop-i
 
1740
c ------- skip g_z1 ---------------------- END
 
1741
c ------- skip g_Az1 ---------------------- START
 
1742
         do ivec=1,nvec
 
1743
          call sread(unitno,dbl_mb(k_zre),n1)
 
1744
          call sread(unitno,dbl_mb(k_zim),n1)
 
1745
         enddo ! end-loop-i
 
1746
c ------- skip g_Az1 ---------------------- END
 
1747
       enddo ! end-loop-iskip
 
1748
c ======= skip blocks ======================= END
 
1749
c ======= Loop in subspaces ===== START
 
1750
       do iblock=imin,imax
 
1751
        read(unitno, err=1001, end=1001) n_read
 
1752
        read(unitno, err=1001, end=1001) ncomp_read
 
1753
        read(unitno, err=1001, end=1001) nvec_read
 
1754
        read(unitno, err=1001, end=1001) n1_read
 
1755
        read(unitno, err=1001, end=1001) nsub_read
 
1756
        write(*,14) iblock,nsub_read
 
1757
 14     format('(iblock,nsub_read)=(',i4,',',i4,')')
 
1758
c ------- read g_z1 ---------------------- START
 
1759
        m1=(iblock-1)*nvec+1
 
1760
        m2=m1+nvec-1
 
1761
c        write(*,4) m1,m2,nvec
 
1762
c 4      format('(m1,m2,nvec)=(',i4,',',i4,',',i4,')')
 
1763
         do ivec=m1,m2 
 
1764
         call dcopy(n1,0.0d0,0,dbl_mb(k_zre),1)
 
1765
         call sread(unitno,dbl_mb(k_zre),n1)
 
1766
         call dcopy(n1,0.0d0,0,dbl_mb(k_zim),1)
 
1767
         call sread(unitno,dbl_mb(k_zim),n1)
 
1768
          do idat=1,n1
 
1769
           val_cmplx=dcmplx(dbl_mb(k_zre+idat-1),
 
1770
     &                      dbl_mb(k_zim+idat-1))
 
1771
           call ga_put(g_z1,idat,idat,ivec,ivec,val_cmplx,1)
 
1772
          enddo ! end-loop-idat
 
1773
         enddo ! end-loop-ivec
 
1774
c ------- read g_z1 ---------------------- END
 
1775
c ------- read g_Az1 ---------------------- START
 
1776
         do ivec=m1,m2
 
1777
          call dcopy(n1,0.0d0,0,dbl_mb(k_zre),1)
 
1778
          call sread(unitno,dbl_mb(k_zre),n1)
 
1779
          call dcopy(n1,0.0d0,0,dbl_mb(k_zim),1)
 
1780
          call sread(unitno,dbl_mb(k_zim),n1)
 
1781
          do idat=1,n1
 
1782
           val_cmplx=dcmplx(dbl_mb(k_zre+idat-1),
 
1783
     &                      dbl_mb(k_zim+idat-1))
 
1784
           call ga_put(g_Az1,idat,idat,ivec,ivec,val_cmplx,1)
 
1785
          enddo ! end-loop-idat
 
1786
         enddo ! end-loop-i
 
1787
c ------- read g_Az1 ---------------------- END
 
1788
        enddo ! end-loop-iblock
 
1789
        nsub=(imax-1)*3
 
1790
        write(*,5) nsub,nblocks_file,nblocks_ga
 
1791
 5      format('dft_CPHF2_read: (nsub,nblocks_file,nblocks_ga)=(',
 
1792
     &         i12,',',i12,',',i12,')')
 
1793
c ======= Loop in subspaces ===== END
 
1794
c     Deallocate the temporary buffer
 
1795
c ----- Using ma_free_heap ------------ START
 
1796
      if (.not.ma_pop_stack(l_zim))
 
1797
     $  call errquit('dft_CPHF2_read: pop problem with l_zim',
 
1798
     &               555,MA_ERR)
 
1799
      if (.not.ma_pop_stack(l_zre))
 
1800
     $  call errquit('dft_CPHF2_read: pop problem with l_zre',
 
1801
     &               555,MA_ERR)
 
1802
c ----- Using ma_free_heap ------------ END
 
1803
c      Close the file
 
1804
       close(unitno,err=1002)
 
1805
       ok = 1
 
1806
      end if ! ----- if-ga_nodeid-eq-0 ----- END
 
1807
c     Broadcast status to other nodes
 
1808
 10   call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status
 
1809
      call ga_sync()
 
1810
      dft_CPHF2_read = ok .eq. 1
 
1811
      return
 
1812
 1000 write(6,*) 'dft_CPHF2_read: failed to open',
 
1813
     $     filename(1:inp_strlen(filename))
 
1814
      call util_flush(luout)
 
1815
      ok = 0
 
1816
      goto 10
 
1817
 1001 write(6,*) 'dft_CPHF2_read: failed to read',
 
1818
     $     filename(1:inp_strlen(filename))
 
1819
      call util_flush(luout)
 
1820
      ok = 0
 
1821
      close(unitno,err=1002)
 
1822
      goto 10
 
1823
 1003 write(6,*) 
 
1824
     & 'dft_CPHF2_read: file inconsistent with calculation',
 
1825
     $     filename(1:inp_strlen(filename))
 
1826
      call util_flush(luout)
 
1827
      ok = 0
 
1828
      close(unitno,err=1002)
 
1829
      goto 10
 
1830
 1002 write(6,*) 'dft_CPHF2_read: failed to close',
 
1831
     $     filename(1:inp_strlen(filename))
 
1832
      call util_flush(luout)
 
1833
      ok = 0
 
1834
      goto 10
 
1835
      end
 
1836
 
 
1837
      logical function dft_CPHF2_read2fix(
 
1838
     &           filename, ! in: filename
 
1839
     &           n,        ! in: sum_{i=1,npol} nocc(i)*nvirt(i)
 
1840
     &           ncomp,    ! in: nr. components
 
1841
     &           nvec,     ! in: nr. of directions = 3
 
1842
     &           n1,       ! in: =n*ncomp
 
1843
     &           nsub,     ! ou: last subspace index (nsub+1)= nr of subspaces stored
 
1844
     &           maxsub,   ! in: max subspace of (g_z1,g_Az1)
 
1845
     &           g_z1,     ! ou: history matrix z
 
1846
     &           g_Az1)    ! ou: history matrix Az
 
1847
c
 
1848
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
 
1849
 
 
1850
      implicit none
 
1851
#include "errquit.fh"
 
1852
#include "global.fh"
 
1853
#include "tcgmsg.fh"
 
1854
#include "msgtypesf.h"
 
1855
#include "mafdecls.fh"
 
1856
#include "msgids.fh"
 
1857
#include "cscfps.fh"
 
1858
#include "inp.fh"
 
1859
#include "util.fh"
 
1860
#include "stdio.fh"
 
1861
      character*(*) filename    ! [input] File to write to
 
1862
      character*255 filename_mini   ! only to store nsub
 
1863
      character*255 filename_mini_1 ! only to store nsub
 
1864
      character*255 filename_1      ! only to store nsub
 
1865
      integer ivec,idat,
 
1866
     &        g_z1,g_Az1
 
1867
      integer unitno
 
1868
      parameter (unitno = 77)
 
1869
 
 
1870
      integer unitno1
 
1871
      parameter (unitno1 = 78)
 
1872
 
 
1873
      integer l_zre,k_zre,
 
1874
     &        l_zim,k_zim
 
1875
      integer ok,i,j,m1,m2,
 
1876
     &        imin,imax,iskip,nskip,nsub1
 
1877
      integer inntsize
 
1878
      integer n,ncomp,nvec,n1,nsub,maxsub,iblock,
 
1879
     &        n_read,n1_read,nvec_red,ncomp_read,
 
1880
     &        nsub_read,nvec_read,
 
1881
     &        nblocks_ga,nblocks_file,
 
1882
     &        ntotblock_true
 
1883
      double complex val_cmplx
 
1884
      external conv2reim_z1 ! defined in ga_lkain_2cpl3.F
 
1885
c 00000000000000000000000000
 
1886
      ntotblock_true=75 ! FePc
 
1887
c 00000000000000000000000000
 
1888
      inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
 
1889
      call ga_sync()
 
1890
      ok = 0
 
1891
      if (ga_nodeid() .eq. 0) then
 
1892
      write(filename_1,31) trim(filename),'_1'
 
1893
 31   format(a,2a)
 
1894
      write(filename_mini_1,32) trim(filename),'_nblock_1'
 
1895
 32   format(a,9a)
 
1896
      write(*,33) filename_1(1:inp_strlen(filename_1)),
 
1897
     &            filename_mini_1(1:inp_strlen(filename_mini_1))
 
1898
 33   format('Creating files: filename_1=',a,
 
1899
     &       ' filename_mini_1=',a)
 
1900
c +++++++++ read nsub +++++++++ START
 
1901
      write(filename_mini,30) trim(filename),'_nblock'
 
1902
 30   format(a,4a)
 
1903
c      write(*,*) 'Reading mini:',filename_mini
 
1904
       open(unitno, status='unknown', form='unformatted',
 
1905
     $      file=filename_mini, err=1000)
 
1906
       read(unitno, err=1001, end=1001) nblocks_file
 
1907
       close(unitno,err=1002)  
 
1908
 
 
1909
      write(*,*) 'Writing nblock=',ntotblock_true
 
1910
       open(unitno1, status='unknown', form='unformatted',
 
1911
     $      file=filename_mini_1, err=1000)
 
1912
       write(unitno1, err=1001) ntotblock_true
 
1913
       close(unitno1,err=1002)  
 
1914
 
 
1915
       nblocks_ga=maxsub/nvec
 
1916
       write(*,2) nblocks_file,nblocks_ga
 
1917
 2     format('(nblocks_file,nblocks_ga)=(',i4,',',i4,')')
 
1918
c --- Compare (nblocks_ga,nblocks_file)
 
1919
       if (nblocks_file .le. nblocks_ga) then
 
1920
        imin=1
 
1921
        imax=nblocks_file
 
1922
        nskip=0
 
1923
       else
 
1924
        imin=1
 
1925
        imax=nblocks_ga
 
1926
        nskip=nblocks_file-nblocks_ga
 
1927
       endif
 
1928
       write(*,3) imin,imax,nskip
 
1929
 3     format('(imin,imax,nskip)=(',i4,',',i4,',',i4,')')
 
1930
c +++++++++ read nsub +++++++++ END
 
1931
c      Print a message indicating the file being read
 
1932
       write(6,22) filename(1:inp_strlen(filename))
 
1933
 22    format(/' dft_CPHF_read:Read aoresponse g_rhs data from ',a/)
 
1934
       call util_flush(luout)
 
1935
c      Open the file
 
1936
       open(unitno, status='old', form='unformatted', file=filename,
 
1937
     $        err=1000)
 
1938
 
 
1939
       open(unitno1, status='unknown', form='unformatted',
 
1940
     $      file=filename_1, err=1000,position='append')
 
1941
 
 
1942
c      Read in some basics to check if they are consistent with the calculation
 
1943
      if (.not.MA_Push_Get(mt_dbl,n1,'hessv jfacs',l_zre,k_zre))
 
1944
     &     call errquit('conv2complex: cannot allocate zre',
 
1945
     &                  n1, MA_ERR)
 
1946
      if (.not.MA_Push_Get(mt_dbl,n1,'hessv kfacs',l_zim,k_zim))
 
1947
     &     call errquit('conv2complex: cannot allocate zim',
 
1948
     &                  n1, MA_ERR)
 
1949
c ======= Loop in subspaces ===== START
 
1950
       nsub1=0
 
1951
       do iblock=1,ntotblock_true
 
1952
        read(unitno, err=1001, end=1001) n_read
 
1953
        read(unitno, err=1001, end=1001) ncomp_read
 
1954
        read(unitno, err=1001, end=1001) nvec_read
 
1955
        read(unitno, err=1001, end=1001) n1_read
 
1956
        read(unitno, err=1001, end=1001) nsub_read
 
1957
 
 
1958
        write(unitno1, err=1001) n_read
 
1959
        write(unitno1, err=1001) ncomp_read
 
1960
        write(unitno1, err=1001) nvec_read
 
1961
        write(unitno1, err=1001) n1_read
 
1962
        write(unitno1, err=1001) nsub1
 
1963
        nsub1=nsub1+nvec
 
1964
        write(*,14) iblock,nsub1,nsub_read
 
1965
 14     format('(iblock,nsub1,nsub_read)=(',
 
1966
     &         i4,',',i4,',',i4,')')
 
1967
c ------- read g_z1 ---------------------- START
 
1968
        m1=(iblock-1)*nvec+1
 
1969
        m2=m1+nvec-1
 
1970
c        write(*,4) m1,m2,nvec
 
1971
c 4      format('(m1,m2,nvec)=(',i4,',',i4,',',i4,')')
 
1972
         do ivec=m1,m2 
 
1973
          call dcopy(n1,0.0d0,0,dbl_mb(k_zre),1)
 
1974
          call sread(unitno,dbl_mb(k_zre),n1)
 
1975
          call dcopy(n1,0.0d0,0,dbl_mb(k_zim),1)
 
1976
          call sread(unitno,dbl_mb(k_zim),n1)
 
1977
          call swrite(unitno1,dbl_mb(k_zre),n1)
 
1978
          call swrite(unitno1,dbl_mb(k_zim),n1)
 
1979
         enddo ! end-loop-i
 
1980
c ------- read g_z1 ---------------------- END
 
1981
c ------- read g_Az1 ---------------------- START
 
1982
         do ivec=m1,m2
 
1983
          call dcopy(n1,0.0d0,0,dbl_mb(k_zre),1)
 
1984
          call sread(unitno,dbl_mb(k_zre),n1)
 
1985
          call dcopy(n1,0.0d0,0,dbl_mb(k_zim),1)
 
1986
          call sread(unitno,dbl_mb(k_zim),n1)
 
1987
          call swrite(unitno1,dbl_mb(k_zre),n1)
 
1988
          call swrite(unitno1,dbl_mb(k_zim),n1)
 
1989
         enddo ! end-loop-i
 
1990
c ------- read g_Az1 ---------------------- END
 
1991
        enddo ! end-loop-iblock
 
1992
c ======= Loop in subspaces ===== END
 
1993
c     Deallocate the temporary buffer
 
1994
c ----- Using ma_free_heap ------------ START
 
1995
      if (.not.ma_pop_stack(l_zim))
 
1996
     $  call errquit('dft_CPHF2_read: pop problem with l_zim',
 
1997
     &               555,MA_ERR)
 
1998
      if (.not.ma_pop_stack(l_zre))
 
1999
     $  call errquit('dft_CPHF2_read: pop problem with l_zre',
 
2000
     &               555,MA_ERR)
 
2001
c ----- Using ma_free_heap ------------ END
 
2002
c      Close the file
 
2003
       close(unitno,err=1002)
 
2004
       close(unitno1,err=1002)
 
2005
       ok = 1
 
2006
      end if
 
2007
c
 
2008
c     Broadcast status to other nodes
 
2009
 10   call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status
 
2010
      call ga_sync()
 
2011
      dft_CPHF2_read2fix = ok .eq. 1
 
2012
      return
 
2013
 1000 write(6,*) 'dft_CPHF2_read: failed to open',
 
2014
     $     filename(1:inp_strlen(filename))
 
2015
      call util_flush(luout)
 
2016
      ok = 0
 
2017
      goto 10
 
2018
 1001 write(6,*) 'dft_CPHF2_read: failed to read',
 
2019
     $     filename(1:inp_strlen(filename))
 
2020
      call util_flush(luout)
 
2021
      ok = 0
 
2022
      close(unitno,err=1002)
 
2023
      goto 10
 
2024
 1003 write(6,*) 
 
2025
     & 'dft_CPHF2_read: file inconsistent with calculation',
 
2026
     $     filename(1:inp_strlen(filename))
 
2027
      call util_flush(luout)
 
2028
      ok = 0
 
2029
      close(unitno,err=1002)
 
2030
      goto 10
 
2031
 1002 write(6,*) 'dft_CPHF2_read: failed to close',
 
2032
     $     filename(1:inp_strlen(filename))
 
2033
      call util_flush(luout)
 
2034
      ok = 0
 
2035
      goto 10
 
2036
      end
 
2037
c =========== READ/WRITE CPHF-1 g_rhs) data ==========END