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
1232
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
1234
c Note.- nmo ne nbf if linear dependencies appear.
1237
#include "errquit.fh"
1238
#include "mafdecls.fh"
1239
#include "global.fh"
1240
#include "tcgmsg.fh"
1241
#include "msgtypesf.h"
1243
#include "msgids.fh"
1244
#include "cscfps.fh"
1247
character*(*) filename ! [input] File to write to
1249
& nocc(npol),nvirt(npol),
1255
parameter (unitno = 77)
1256
integer l_rhs0,k_rhs0,
1264
inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
1267
c Read routines should be consistent with this
1268
c Write out the atomic zora corrections
1269
if (ga_nodeid() .eq. 0) then
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
1276
write(unitno, err=1001) nocc(i)
1279
write(unitno, err=1001) nvirt(i)
1281
write(unitno, err=1001) nxyz
1282
c Allocate the temporary buffer
1285
ntot=ntot+nocc(ispin)*nvirt(ispin)
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',
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)
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)
1304
endif ! end-if-lifetime
1305
enddo ! end-loop-ipm
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',
1312
c ----- Using ma_free_heap ------------ END
1314
close(unitno,err=1002)
1317
c Broadcast status to other nodes
1318
10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status
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)
1327
1000 write(6,*) 'dft_CPHF_write: failed to open ',
1328
$ filename(1:inp_strlen(filename))
1329
call util_flush(luout)
1332
1001 write(6,*) 'dft_CPHF_write: failed to write ',
1333
$ filename(1:inp_strlen(filename))
1334
call util_flush(luout)
1336
close(unitno,err=1002)
1338
1002 write(6,*) 'dft_CPHF_write: failed to close',
1339
$ filename(1:inp_strlen(filename))
1340
call util_flush(luout)
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
1355
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
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"
1368
character*(*) filename ! [input] File to write to
1370
& nocc(npol),nvirt(npol),
1376
parameter (unitno = 77)
1381
integer npol_read,nxyz_read,ntot_read,
1382
& nocc_read(2),nvirt_read(2)
1385
c Initialise to invalid MA handle
1386
inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
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)
1395
open(unitno, status='old', form='unformatted', file=filename,
1397
c Read in some basics to check if they are consistent with the calculation
1398
read(unitno, err=1001, end=1001) npol_read
1400
read(unitno, err=1001, end=1001) nocc_read(i)
1403
read(unitno, err=1001, end=1001) nvirt_read(i)
1405
read(unitno, err=1001, end=1001) nxyz_read
1407
if ((nxyz_read .ne. nxyz) .or.
1408
& (npol_read .ne. npol)) goto 1003
1411
ntot=ntot+nocc(ispin)*nvirt(ispin)
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',
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)
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)
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)
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)
1439
endif ! end-if-lifetime
1440
enddo ! end-loop-ipm
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',
1447
c ----- Using ma_free_heap ------------ END
1449
close(unitno,err=1002)
1453
c Broadcast status to other nodes
1454
10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status
1456
dft_CPHF1_read = ok .eq. 1
1458
1000 write(6,*) 'dft_CPHF_read: failed to open',
1459
$ filename(1:inp_strlen(filename))
1460
call util_flush(luout)
1463
1001 write(6,*) 'dft_CPHF_read: failed to read',
1464
$ filename(1:inp_strlen(filename))
1465
call util_flush(luout)
1467
close(unitno,err=1002)
1470
& 'dft_CPHF_read: file inconsistent with calculation',
1471
$ filename(1:inp_strlen(filename))
1472
call util_flush(luout)
1474
close(unitno,err=1002)
1476
1002 write(6,*) 'dft_CPHF_read: failed to close',
1477
$ filename(1:inp_strlen(filename))
1478
call util_flush(luout)
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
1495
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
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
1501
c -->It will write only subscape nsub for (z1,Az1)
1503
#include "errquit.fh"
1504
#include "mafdecls.fh"
1505
#include "global.fh"
1506
#include "tcgmsg.fh"
1507
#include "msgtypesf.h"
1509
#include "msgids.fh"
1510
#include "cscfps.fh"
1513
character*(*) filename ! [input] File to write to
1514
character*255 filename_mini ! only to store nblocks
1517
parameter (unitno = 77)
1518
integer l_dat,k_dat,
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
1526
inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
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'
1535
c write(*,*) 'Writing mini:',filename_mini
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
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',
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',
1563
c ------- write g_z1 ---------------------- START
1565
m1=(nblock-1)*nvec+i
1566
call ga_get(g_z1,1,n1,m1,m1,dcpl_mb(k_z),1)
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)
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)
1578
c ------- write g_z1 ---------------------- END
1579
c ------- write g_Az1 ---------------------- START
1581
m1=(nblock-1)*nvec+i
1582
call ga_get(g_Az1,1,n1,m1,m1,dcpl_mb(k_z),1)
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)
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)
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',
1600
if (.not. ma_free_heap(l_z))
1601
$ call errquit('dft_CPHF2_write: ma free_heap failed',
1603
c ----- Using ma_free_heap ------------ END
1605
close(unitno,err=1002)
1608
c Broadcast status to other nodes
1609
10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status
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)
1622
1000 write(6,*) 'dft_CPHF2_write: failed to open ',
1623
$ filename(1:inp_strlen(filename))
1624
call util_flush(luout)
1627
1001 write(6,*) 'dft_CPHF2_write: failed to write ',
1628
$ filename(1:inp_strlen(filename))
1629
call util_flush(luout)
1631
close(unitno,err=1002)
1633
1002 write(6,*) 'dft_CPHF2_write: failed to close',
1634
$ filename(1:inp_strlen(filename))
1635
call util_flush(luout)
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
1652
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
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"
1665
character*(*) filename ! [input] File to write to
1666
character*255 filename_mini ! only to store nsub
1670
parameter (unitno = 77)
1671
integer l_zre,k_zre,
1673
integer ok,i,j,m1,m2,
1674
& imin,imax,iskip,nskip
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
1684
inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
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'
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
1709
nskip=nblocks_file-(nblocks_ga-1)
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)
1719
open(unitno, status='old', form='unformatted', file=filename,
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',
1725
if (.not.MA_Push_Get(mt_dbl,n1,'hessv kfacs',l_zim,k_zim))
1726
& call errquit('conv2complex: cannot allocate zim',
1728
c ======= skip blocks ======================= START
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
1737
call sread(unitno,dbl_mb(k_zre),n1)
1738
call sread(unitno,dbl_mb(k_zim),n1)
1740
c ------- skip g_z1 ---------------------- END
1741
c ------- skip g_Az1 ---------------------- START
1743
call sread(unitno,dbl_mb(k_zre),n1)
1744
call sread(unitno,dbl_mb(k_zim),n1)
1746
c ------- skip g_Az1 ---------------------- END
1747
enddo ! end-loop-iskip
1748
c ======= skip blocks ======================= END
1749
c ======= Loop in subspaces ===== START
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
1761
c write(*,4) m1,m2,nvec
1762
c 4 format('(m1,m2,nvec)=(',i4,',',i4,',',i4,')')
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)
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
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)
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
1787
c ------- read g_Az1 ---------------------- END
1788
enddo ! end-loop-iblock
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',
1799
if (.not.ma_pop_stack(l_zre))
1800
$ call errquit('dft_CPHF2_read: pop problem with l_zre',
1802
c ----- Using ma_free_heap ------------ END
1804
close(unitno,err=1002)
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
1810
dft_CPHF2_read = ok .eq. 1
1812
1000 write(6,*) 'dft_CPHF2_read: failed to open',
1813
$ filename(1:inp_strlen(filename))
1814
call util_flush(luout)
1817
1001 write(6,*) 'dft_CPHF2_read: failed to read',
1818
$ filename(1:inp_strlen(filename))
1819
call util_flush(luout)
1821
close(unitno,err=1002)
1824
& 'dft_CPHF2_read: file inconsistent with calculation',
1825
$ filename(1:inp_strlen(filename))
1826
call util_flush(luout)
1828
close(unitno,err=1002)
1830
1002 write(6,*) 'dft_CPHF2_read: failed to close',
1831
$ filename(1:inp_strlen(filename))
1832
call util_flush(luout)
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
1848
c Author : Fredy W. Aquino, Northwestern University (Oct 2012)
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"
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
1868
parameter (unitno = 77)
1871
parameter (unitno1 = 78)
1873
integer l_zre,k_zre,
1875
integer ok,i,j,m1,m2,
1876
& imin,imax,iskip,nskip,nsub1
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,
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)
1891
if (ga_nodeid() .eq. 0) then
1892
write(filename_1,31) trim(filename),'_1'
1894
write(filename_mini_1,32) trim(filename),'_nblock_1'
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'
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)
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)
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
1926
nskip=nblocks_file-nblocks_ga
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)
1936
open(unitno, status='old', form='unformatted', file=filename,
1939
open(unitno1, status='unknown', form='unformatted',
1940
$ file=filename_1, err=1000,position='append')
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',
1946
if (.not.MA_Push_Get(mt_dbl,n1,'hessv kfacs',l_zim,k_zim))
1947
& call errquit('conv2complex: cannot allocate zim',
1949
c ======= Loop in subspaces ===== START
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
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
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
1970
c write(*,4) m1,m2,nvec
1971
c 4 format('(m1,m2,nvec)=(',i4,',',i4,',',i4,')')
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)
1980
c ------- read g_z1 ---------------------- END
1981
c ------- read g_Az1 ---------------------- START
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)
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',
1998
if (.not.ma_pop_stack(l_zre))
1999
$ call errquit('dft_CPHF2_read: pop problem with l_zre',
2001
c ----- Using ma_free_heap ------------ END
2003
close(unitno,err=1002)
2004
close(unitno1,err=1002)
2008
c Broadcast status to other nodes
2009
10 call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status
2011
dft_CPHF2_read2fix = ok .eq. 1
2013
1000 write(6,*) 'dft_CPHF2_read: failed to open',
2014
$ filename(1:inp_strlen(filename))
2015
call util_flush(luout)
2018
1001 write(6,*) 'dft_CPHF2_read: failed to read',
2019
$ filename(1:inp_strlen(filename))
2020
call util_flush(luout)
2022
close(unitno,err=1002)
2025
& 'dft_CPHF2_read: file inconsistent with calculation',
2026
$ filename(1:inp_strlen(filename))
2027
call util_flush(luout)
2029
close(unitno,err=1002)
2031
1002 write(6,*) 'dft_CPHF2_read: failed to close',
2032
$ filename(1:inp_strlen(filename))
2033
call util_flush(luout)
2037
c =========== READ/WRITE CPHF-1 g_rhs) data ==========END