← Back to branch summary

~ubuntu-branches/ubuntu/trusty/nwchem/trusty-proposed

« back to all changes in this revision

Viewing changes to src/nwpw/pspw/lib/psp/psp.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/nwpw/pspw/lib/psp/psp.F
1
1
*
2
 
* $Id: psp.F 21861 2012-01-25 20:08:42Z bylaska $
 
2
* $Id: psp.F 23826 2013-03-18 18:04:02Z bylaska $
3
3
*
4
4
 
5
5
*     ***********************************
22
22
      integer taskid,MASTER
23
23
      parameter (MASTER=0)
24
24
 
25
 
      integer npack1,npack0,nion,i
 
25
      integer npack1,npack0,nion,i,j,ii
26
26
      logical value
27
27
 
28
28
*     **** external functions *****
29
 
      logical  control_pspspin
30
 
      external control_pspspin
31
 
      integer  ion_nkatm,ion_nion
32
 
      external ion_nkatm,ion_nion
 
29
      logical  control_pspspin,control_use_grid_cmp
 
30
      external control_pspspin,control_use_grid_cmp
 
31
      integer  ion_nkatm,ion_nion,control_pspnuterms
 
32
      external ion_nkatm,ion_nion,control_pspnuterms
33
33
 
34
34
      call Pack_npack(1,npack1)
35
35
      call Pack_npack(0,npack0)
45
45
c    >                    'vnlnrm',vnlnrm(2),vnlnrm(1))
46
46
 
47
47
      value = MA_alloc_get(mt_dbl,(npsp*npack0),'vl',vl(2),vl(1))
48
 
      value = value.and.MA_alloc_get(mt_int,npsp,'vnl',vnl(2),vnl(1))
 
48
      value = value.and.
 
49
     >        MA_alloc_get(mt_int,npsp,'vnl',vnl(2),vnl(1))
49
50
c      value = value.and.
50
51
c     >      MA_alloc_get(mt_dbl,(nmax_max*nmax_max*(lmax_max+1)*npsp),
51
52
c     >                    'Gijl',Gijl(2),Gijl(1))
186
187
*     **** allocate semicore data ****
187
188
      call semicore_init()
188
189
 
189
 
*     *** set pawexist ****
 
190
*     *** set pawexist and use_grid_cmp ****
190
191
      pawexist = .false.
 
192
      use_grid_cmp = control_use_grid_cmp()
 
193
 
191
194
 
192
195
      return
193
196
      end
206
209
 
207
210
*     ***** local variables ****
208
211
      logical value
209
 
      integer npack1
 
212
      integer npack1,npack0
210
213
 
211
214
*     ***** external functions *****
212
215
      integer  psp_nprj_max,control_nprj_mult
213
216
      external psp_nprj_max,control_nprj_mult
214
217
 
215
218
 
 
219
      call Pack_npack(0,npack0)
216
220
      call Pack_npack(1,npack1)
217
221
      nprj_mult= control_nprj_mult()
218
222
      nprj_max = psp_nprj_max()
219
223
      value = MA_alloc_get(mt_dcpl,npack1*nprj_max*nprj_mult,
220
224
     >                     'prjtmp',prjtmp(2),prjtmp(1))
 
225
      if (pawexist) then
 
226
         value = value.and.
 
227
     >            MA_alloc_get(mt_dbl,2*nprj_max*nprj_max,
 
228
     >                        'wtmp',wtmp(2),wtmp(1))
 
229
         value = value.and.
 
230
     >            MA_alloc_get(mt_dcpl,npack0,'vc_tmp',
 
231
     >                         vc_tmp(2),vc_tmp(1))
 
232
         value = value.and.
 
233
     >            MA_alloc_get(mt_dcpl,npack0,'vcmptmp',
 
234
     >                         vcmp_tmp(2),vcmp_tmp(1))
 
235
         value = value.and.
 
236
     >            MA_alloc_get(mt_dcpl,npack0,'vlpaw_tmp',
 
237
     >                         vlpaw_tmp(2),vlpaw_tmp(1))
 
238
      end if
221
239
      if (.not. value)
222
240
     >   call errquit('psp_proj_init:out of heap memory',0,MA_ERR)
223
241
      return
224
242
      end
225
243
 
 
244
      subroutine psp_set_vcandvl(vcin,vlin)
 
245
      implicit none
 
246
      complex*16 vcin(*)
 
247
      complex*16 vlin(*)
 
248
 
 
249
#include "errquit.fh"
 
250
#include "mafdecls.fh"
 
251
#include "psp.fh"
 
252
 
 
253
*     **** local variables ****
 
254
      real*8 scal
 
255
 
 
256
*     **** external functions ****
 
257
      real*8   lattice_omega
 
258
      external lattice_omega
 
259
 
 
260
      scal = 1.0d0/lattice_omega()
 
261
 
 
262
      call Pack_c_Copy(0,vcin,dcpl_mb(vc_tmp(1))) 
 
263
      call Pack_cc_daxpy(0,scal,vlin,dcpl_mb(vc_tmp(1))) 
 
264
      return
 
265
      end 
 
266
 
 
267
      subroutine psp_get_vcmpvlpaw(vcmpout)
 
268
      implicit none
 
269
      complex*16 vcmpout(*)
 
270
 
 
271
#include "errquit.fh"
 
272
#include "mafdecls.fh"
 
273
#include "psp.fh"
 
274
 
 
275
*     **** local variables ****
 
276
      real*8 scal
 
277
 
 
278
*     **** external functions ****
 
279
      real*8   lattice_omega
 
280
      external lattice_omega
 
281
 
 
282
      scal = 1.0d0/lattice_omega()
 
283
 
 
284
      call Pack_c_Copy(0,dcpl_mb(vcmp_tmp(1)),vcmpout) 
 
285
      call Pack_cc_daxpy(0,scal,dcpl_mb(vlpaw_tmp(1)),vcmpout) 
 
286
      return
 
287
      end 
 
288
 
 
289
 
 
290
 
 
291
      real*8 function psp_vlpaw_e(dng)
 
292
      implicit none
 
293
      complex*16 dng(*)
 
294
 
 
295
#include "errquit.fh"
 
296
#include "mafdecls.fh"
 
297
#include "psp.fh"
 
298
 
 
299
*     **** local variables ****
 
300
      real*8 e
 
301
 
 
302
*     **** external functions ****
 
303
      real*8   lattice_omega
 
304
      external lattice_omega
 
305
 
 
306
      call Pack_cc_dot(0,dcpl_mb(vlpaw_tmp(1)),dng,e)
 
307
      psp_vlpaw_e = e
 
308
      return
 
309
      end
 
310
 
 
311
 
 
312
 
226
313
 
227
314
*     ***********************************
228
315
*     *                                 *
277
364
      end do
278
365
 
279
366
      value = MA_free_heap(prjtmp(2))
 
367
      if (pawexist) then
 
368
         value = value.and.MA_free_heap(wtmp(2))
 
369
         value = value.and.MA_free_heap(vc_tmp(2))
 
370
         value = value.and.MA_free_heap(vcmp_tmp(2))
 
371
         value = value.and.MA_free_heap(vlpaw_tmp(2))
 
372
      end if
280
373
      value = value.and.MA_free_heap(vl(2))
281
374
      value = value.and.MA_free_heap(vnl(2))
282
375
c     value = value.and.MA_free_heap(vnlnrm(2))
334
427
         value = value.and.MA_free_heap(pspspin_downions(2))
335
428
      end if
336
429
 
 
430
 
337
431
      if (.not. value) 
338
432
     >  call errquit('psp_end:error freeing heap memory',0,MA_ERR)
339
433
 
340
 
      return
341
 
      end
 
434
*     **** deallocate prj_indx ****
 
435
      call psp_prj_indx_end()
 
436
 
 
437
      return
 
438
      end
 
439
 
 
440
 
 
441
*     ****************************************************
 
442
*     *                                                  *
 
443
*     *             psp_prj_indx_init                    *
 
444
*     *                                                  *
 
445
*     ****************************************************
 
446
*
 
447
*     This routine sets up the prj_indx indexes:
 
448
*           shift_prj_indx,ii_prj_indx,ia_prj_indx,sd_function_prj_indx
 
449
*
 
450
      subroutine psp_prj_indx_init()
 
451
      implicit none
 
452
 
 
453
#include "mafdecls.fh"
 
454
#include "psp.fh"
 
455
#include "errquit.fh"
 
456
 
 
457
*     **** local variables ****
 
458
      logical value,sd_function
 
459
      integer k,l,ii,ia,nproj,l_prj,count1,count2,shift,nion
 
460
 
 
461
*     **** external functions ****
 
462
      integer  ion_nion,ion_katm,psi_data_get_ptr
 
463
      external ion_nion,ion_katm,psi_data_get_ptr
 
464
 
 
465
      nion = ion_nion()
 
466
      n_prj_indx    = 0
 
467
      nion_prj_indx = 0
 
468
      do ii=1,nion
 
469
        ia=ion_katm(ii)
 
470
 
 
471
        nproj = int_mb(nprj(1)+ia-1)
 
472
 
 
473
        if (nproj.gt.0) then
 
474
           nion_prj_indx = nion_prj_indx + 1
 
475
 
 
476
        do l=1,nproj
 
477
           n_prj_indx = n_prj_indx + 1
 
478
        end do
 
479
        end if !** nproj>0 **
 
480
      end do !** ii **
 
481
 
 
482
      value = MA_alloc_get(mt_int,nion_prj_indx,'ii_prj_indx',
 
483
     >                     ii_prj_indx(2),ii_prj_indx(1))
 
484
      value = value.and.
 
485
     >        MA_alloc_get(mt_int,nion_prj_indx,'ia_prj_indx',
 
486
     >                     ia_prj_indx(2),ia_prj_indx(1))
 
487
      value = value.and.
 
488
     >        MA_alloc_get(mt_int,nion_prj_indx,'nproj_prj_indx',
 
489
     >                     nproj_prj_indx(2),nproj_prj_indx(1))
 
490
      value = value.and.
 
491
     >        MA_alloc_get(mt_int,nion_prj_indx,'swstart_prj_indx',
 
492
     >                     swstart_prj_indx(2),swstart_prj_indx(1))
 
493
 
 
494
      value = MA_alloc_get(mt_int,n_prj_indx,'shift_prj_indx',
 
495
     >                     shift_prj_indx(2),shift_prj_indx(1))
 
496
      value = MA_alloc_get(mt_int,n_prj_indx,'l_prj_prj_indx',
 
497
     >                     l_prj_prj_indx(2),l_prj_prj_indx(1))
 
498
      value = value.and.
 
499
     >        MA_alloc_get(mt_log,n_prj_indx,'sd_function_prj_indx',
 
500
     >              sd_function_prj_indx(2),sd_function_prj_indx(1))
 
501
      if (.not.value)
 
502
     >  call errquit('psp_prj_indx_init: out of memory',0, MA_ERR)
 
503
      swaset = .false.
 
504
      swann  = 0
 
505
 
 
506
      count1 = 0
 
507
      count2 = 0
 
508
      do ii=1,nion
 
509
        ia    = ion_katm(ii)
 
510
        nproj = int_mb(nprj(1)+ia-1)
 
511
        if (nproj.gt.0) then
 
512
           int_mb(ii_prj_indx(1)+count1)      = ii
 
513
           int_mb(ia_prj_indx(1)+count1)      = ia
 
514
           int_mb(nproj_prj_indx(1)+count1)   = nproj
 
515
           int_mb(swstart_prj_indx(1)+count1) = count2
 
516
           count1 = count1 + 1
 
517
 
 
518
           do l=1,nproj
 
519
              shift = psi_data_get_ptr(int_mb(vnl(1)+ia-1),l)
 
520
              l_prj = int_mb(l_projector(1)+(l-1)
 
521
     >                                  + (ia-1)*(nmax_max*lmmax_max))
 
522
#ifdef GCC4
 
523
              k = iand(l_prj,1)
 
524
#else
 
525
              k = and(l_prj,1)
 
526
#endif
 
527
              sd_function = (k.eq.0)
 
528
 
 
529
              int_mb(shift_prj_indx(1)+count2)       = shift
 
530
              int_mb(l_prj_prj_indx(1)+count2)       = l_prj
 
531
              log_mb(sd_function_prj_indx(1)+count2) = sd_function
 
532
 
 
533
              count2 = count2 + 1
 
534
           end do
 
535
        end if !** nproj>0 **
 
536
      end do !** ii **
 
537
 
 
538
      return 
 
539
      end
 
540
 
 
541
 
 
542
*     ****************************************************
 
543
*     *                                                  *
 
544
*     *         psp_prj_indx_alloc_sw1a_sw2a             *
 
545
*     *                                                  *
 
546
*     ****************************************************
 
547
      subroutine psp_prj_indx_alloc_sw1a_sw2a(nn)
 
548
      implicit none
 
549
      integer nn
 
550
 
 
551
#include "mafdecls.fh"
 
552
#include "psp.fh"
 
553
#include "errquit.fh"
 
554
 
 
555
*     **** local variables ****
 
556
      logical ok
 
557
 
 
558
      if (swaset.and.(nn.gt.swann)) then
 
559
         ok =        MA_free_heap(sw1a(2))
 
560
         ok = ok.and.MA_free_heap(sw2a(2))
 
561
         if (.not.ok)
 
562
     >     call errquit('psp_prj_indx_alloc_sw1a_sw2a:freeing heap',
 
563
     >                  0,MA_ERR)
 
564
         swaset = .false.
 
565
      end if
 
566
 
 
567
      if (.not.swaset) then
 
568
         ok = MA_alloc_get(mt_dbl,nn*n_prj_indx,'sw1a',sw1a(2),sw1a(1))
 
569
         ok = ok.and.
 
570
     >        MA_alloc_get(mt_dbl,nn*n_prj_indx,'sw2a',sw2a(2),sw2a(1))
 
571
         if (.not.ok)
 
572
     >     call errquit('psp_prj_indx_alloc_sw1a_sw2a: out of memory',
 
573
     >                  0,MA_ERR)
 
574
         swaset = .true.
 
575
         swann  = nn
 
576
      end if
 
577
 
 
578
      return
 
579
      end
 
580
 
 
581
 
 
582
 
 
583
*     ****************************************************
 
584
*     *                                                  *
 
585
*     *             psp_prj_indx_end                     *
 
586
*     *                                                  *
 
587
*     ****************************************************
 
588
      subroutine psp_prj_indx_end()
 
589
      implicit none
 
590
 
 
591
#include "mafdecls.fh"
 
592
#include "psp.fh"
 
593
#include "errquit.fh"
 
594
 
 
595
*     **** local variables ****
 
596
      logical value
 
597
 
 
598
      value =           MA_free_heap(ii_prj_indx(2))
 
599
      value = value.and.MA_free_heap(ia_prj_indx(2))
 
600
      value = value.and.MA_free_heap(nproj_prj_indx(2))
 
601
      value = value.and.MA_free_heap(swstart_prj_indx(2))
 
602
      value = value.and.MA_free_heap(shift_prj_indx(2))
 
603
      value = value.and.MA_free_heap(l_prj_prj_indx(2))
 
604
      value = value.and.MA_free_heap(sd_function_prj_indx(2))
 
605
      if (swaset) then
 
606
         value = value.and.MA_free_heap(sw1a(2))
 
607
         value = value.and.MA_free_heap(sw2a(2))
 
608
      end if
 
609
      if (.not.value)
 
610
     >  call errquit('psp_prj_indx_end:error freeing heap',0,MA_ERR)
 
611
 
 
612
      return
 
613
      end
 
614
 
342
615
 
343
616
 
344
617
*     ***********************************
809
1082
      end
810
1083
 
811
1084
 
812
 
 
813
1085
*     **********************************
814
1086
*     *                                *
815
1087
*     *        grad_v_lr_local         *
1074
1346
         call dcopy(3*nion,0.0d0,0,fion,1)
1075
1347
      end if
1076
1348
 
 
1349
      if (pawexist) 
 
1350
     >   call dcopy((2*npack0),0.0d0,0,dcpl_mb(vlpaw_tmp(1)),1)
 
1351
 
1077
1352
      call dcopy((2*npack0),0.0d0,0,vl_out,1)
1078
1353
      do ii=1,nion
1079
1354
       if (mod(ii-1,np_j).eq.taskid_j) then
1100
1375
           call Pack_tc_Mul(0,dbl_mb(vl(1)+npack0*(ia-1)),
1101
1376
     >                      dcpl_mb(exi(1)),
1102
1377
     >                      dcpl_mb(vtmp(1)))
1103
 
           call Pack_cc_Sum2(0,dcpl_mb(vtmp(1)),vl_out)
 
1378
 
 
1379
           if ((int_mb(psp_type(1)+ia-1).eq.4)) then
 
1380
              call Pack_cc_Sum2(0,dcpl_mb(vtmp(1)),
 
1381
     >                            dcpl_mb(vlpaw_tmp(1)))
 
1382
           else
 
1383
              call Pack_cc_Sum2(0,dcpl_mb(vtmp(1)),vl_out)
 
1384
           end if
1104
1385
 
1105
1386
           if (move) then
1106
1387
#ifndef CRAY
1119
1400
 
1120
1401
       end if
1121
1402
      end do
1122
 
      if (np_j.gt.1) call D1dB_Vector_SumAll(2*npack0,vl_out)
 
1403
      if (np_j.gt.1) then
 
1404
         call D1dB_Vector_SumAll(2*npack0,vl_out)
 
1405
         if (pawexist) 
 
1406
     >      call D1dB_Vector_SumAll(2*npack0,dcpl_mb(vlpaw_tmp(2)))
 
1407
      end if
1123
1408
      if (move) call Parallel_Vector_SumAll(3*nion,fion)
1124
1409
 
1125
1410
      value = .true.
1239
1524
 
1240
1525
*     ***********************************
1241
1526
*     *                                 *
1242
 
*     *            v_nonlocal           *
 
1527
*     *            v_nonlocal_old       *
1243
1528
*     *                                 *
1244
1529
*     ***********************************
1245
1530
 
1249
1534
*  Note - This routine was restructured 5-13-2002 to improve
1250
1535
*         parallel efficiency.
1251
1536
*
1252
 
      subroutine v_nonlocal(ispin,ne,psi1,psi2,move,fion,fractional,occ)
 
1537
      subroutine v_nonlocal_old(ispin,ne,psi1,psi2,move,fion,
 
1538
     >                          fractional,occ)
1253
1539
      implicit none
1254
1540
      integer    ispin,ne(2)
1255
1541
      complex*16 psi1(*)
1262
1548
#include "mafdecls.fh"
1263
1549
#include "psp.fh"
1264
1550
#include "errquit.fh"
1265
 
ccccccc#include "frac_occ.fh"
1266
1551
 
1267
1552
 
1268
1553
*     *** local variables ***
1269
 
      integer G(3),npack1,nion
1270
 
      integer i,ii,ia,l,n,nn
 
1554
      integer G(3),npack1,nion,nu
 
1555
      integer i,j,ii,ia,l,n,nn
1271
1556
      integer k,shift,l_prj,nproj,Gijl_indx
1272
1557
      real*8  omega,scal,ff(3)
1273
1558
      complex*16 ctmp
1274
 
      integer exi(2),xtmp(2),sw1(2),sw2(2),sum(2)
 
1559
      integer exi(2),xtmp(2),sw1(2),sw2(2),sw3(2),sum(2)
1275
1560
      logical value,sd_function
 
1561
      real*8 vmm(50)
 
1562
      integer ld_ptr
1276
1563
 
1277
1564
*     **** external functions ****
1278
1565
      logical  is_sORd
1380
1667
     >                       dbl_mb(sw1(1)+(l-1)*nn+ne(1)),1)
1381
1668
            end if
1382
1669
        end do
1383
 
        call D3dB_Vector_SumAll((nn*nproj*nprj_mult),dbl_mb(sw1(1)))
 
1670
        call D3dB_Vector_SumAll((nn*nproj),dbl_mb(sw1(1)))
1384
1671
 
1385
1672
 
1386
1673
*       **** sw2 = Gijl*sw1 ******
1398
1685
     >                         dbl_mb(Gijl_indx),
1399
1686
     >                         dbl_mb(sw1(1)),
1400
1687
     >                         dbl_mb(sw2(1)))
1401
 
        
 
1688
 
1402
1689
*       **** do Kleinman-Bylander Multiplication ****
1403
1690
        !scal = 1.0d0/(omega)
1404
 
        call dscal(nn*int_mb(nprj(1)+ia-1)*nprj_mult,
 
1691
        call dscal(nn*int_mb(nprj(1)+ia-1),
1405
1692
     >             scal,dbl_mb(sw2(1)),1)
1406
1693
 
1407
1694
*       **** add xc and coulomb paw parts to sw2 ***
1426
1713
 
1427
1714
        end if
1428
1715
 
1429
 
        call DGEMM('N','T',2*npack1,nn,int_mb(nprj(1)+ia-1)*nprj_mult,
 
1716
        call DGEMM('N','T',2*npack1,nn,int_mb(nprj(1)+ia-1),
1430
1717
     >             (-1.0d0),
1431
1718
     >             dcpl_mb(prjtmp(1)), 2*npack1,
1432
1719
     >             dbl_mb(sw2(1)),     nn,
1504
1791
      if (.not.value) call errquit('v_nonlocal: popping stack',3,
1505
1792
     &       MA_ERR)
1506
1793
      call nwpw_timing_end(6)
1507
 
      return 
1508
 
      end
 
1794
 
 
1795
      return 
 
1796
      end
 
1797
 
 
1798
 
 
1799
*     ***********************************
 
1800
*     *                                 *
 
1801
*     *        v_nonlocal               *
 
1802
*     *                                 *
 
1803
*     ***********************************
 
1804
 
 
1805
*    This routine computes the Kleinman-Bylander non-local 
 
1806
* pseudopotential projection.
 
1807
*
 
1808
*  Note - This routine was restructured 12-1-2013 to handle PAW operators.
 
1809
*
 
1810
*  To Do -  For very large numbers of atoms the code will need to be restructured
 
1811
*           to distribute the sw1a and sw2a matrices over np_i.  Basically, if the orthogonalization matrices
 
1812
*           need to be distributed then sw1a and sw2a will need to be distributed as well.  A simple algorithm
 
1813
*           to do this will be to keep the loop structure the same but instead of the
 
1814
*           call to D3dB_Vector_SumAll(nn*n_prj_indx,dbl_mb(sw1a(1))), this will need to be, 
 
1815
*           changed to D3dB_Vector_SumAll(nn*nproj,dbl_mb(sw1(1))) and then place sw1 --> sw1a(distributed)
 
1816
*
 
1817
      subroutine v_nonlocal(ispin,ne,psi1,psi2,move,fion,
 
1818
     >                      fractional,occ)
 
1819
      implicit none
 
1820
      integer    ispin,ne(2)
 
1821
      complex*16 psi1(*)
 
1822
      complex*16 psi2(*)
 
1823
      logical move
 
1824
      real*8 fion(3,*)
 
1825
      logical fractional
 
1826
      real*8 occ(*)
 
1827
 
 
1828
#include "mafdecls.fh"
 
1829
#include "psp.fh"
 
1830
#include "errquit.fh"
 
1831
 
 
1832
*     *** local variables ***
 
1833
      integer pcount,taskid_j,np_j
 
1834
      integer G(3),npack0,npack1,nion,nu,mult_l,m,ms
 
1835
      integer i,j,ii,ia,l,n,nn,l1,ip,iip,iipmax,jp,swstart
 
1836
      integer k,shift,l_prj,nproj,Gijl_indx
 
1837
      integer nx,ny,nz
 
1838
      real*8  omega,scal,ff(3),scal1
 
1839
      complex*16 ctmp
 
1840
      integer exi(2),xtmp(2),sum(2)
 
1841
      integer dng_cmp(2),dng_cmp_smooth(2)
 
1842
      integer vcmp(2),vcmp_smooth(2)
 
1843
      logical value,sd_function,periodic
 
1844
      real*8 vmm(50),eh_atom
 
1845
      integer ld_ptr
 
1846
 
 
1847
*     **** external functions ****
 
1848
      logical  is_sORd
 
1849
      integer  ion_nion,ion_katm,Pack_G_indx
 
1850
      integer  psi_data_get_ptr,psi_data_get_chnk
 
1851
      real*8   lattice_omega,ddot
 
1852
      external is_sORd
 
1853
      external ion_nion,ion_katm,Pack_G_indx
 
1854
      external psi_data_get_ptr,psi_data_get_chnk
 
1855
      external lattice_omega,ddot
 
1856
      integer  nwpw_compcharge_mult_l,control_version
 
1857
      external nwpw_compcharge_mult_l,control_version
 
1858
      real*8   nwpw_compcharge_Qlm
 
1859
      external nwpw_compcharge_Qlm
 
1860
 
 
1861
      call nwpw_timing_start(6) 
 
1862
 
 
1863
      periodic = (control_version().eq.3)
 
1864
      call Parallel2d_taskid_j(taskid_j)
 
1865
      call Parallel2d_np_j(np_j)
 
1866
 
 
1867
*     **** allocate local memory ****
 
1868
      nion = ion_nion()
 
1869
      nn = ne(1)+ne(2)
 
1870
      call Pack_npack(1,npack1)
 
1871
 
 
1872
      call psp_prj_indx_alloc_sw1a_sw2a(nn)
 
1873
      value = MA_push_get(mt_dcpl,npack1,'exi',exi(2), exi(1))
 
1874
      if (.not.value) 
 
1875
     >  call errquit('v_nonlocal:out of stack',0, MA_ERR)
 
1876
 
 
1877
      if (move) then
 
1878
       value = value.and.MA_push_get(mt_dbl,npack1,
 
1879
     >                               'xtmp',xtmp(2),xtmp(1))
 
1880
       value = value.and.MA_push_get(mt_dbl,3*nn,'sum',sum(2),sum(1))
 
1881
       if (.not. value) 
 
1882
     >  call errquit('v_nonlocal:out of stack',1,MA_ERR)
 
1883
 
 
1884
       G(1)  = Pack_G_indx(1,1)
 
1885
       G(2)  = Pack_G_indx(1,2)
 
1886
       G(3)  = Pack_G_indx(1,3)
 
1887
      end if
 
1888
 
 
1889
      omega = lattice_omega()
 
1890
      scal = 1.0d0/(omega)
 
1891
 
 
1892
 
 
1893
      jp = 0
 
1894
      do ip=1,nion_prj_indx
 
1895
         ii          = int_mb(ii_prj_indx(1)+ip-1)
 
1896
         ia          = int_mb(ia_prj_indx(1)+ip-1)
 
1897
         nproj       = int_mb(nproj_prj_indx(1)+ip-1)
 
1898
 
 
1899
*        **** structure factor and local pseudopotential ****
 
1900
         call strfac_pack(1,ii,dcpl_mb(exi(1)))
 
1901
 
 
1902
         do l=1,nproj
 
1903
            shift       = int_mb(shift_prj_indx(1)+jp)
 
1904
            sd_function = log_mb(sd_function_prj_indx(1)+jp)
 
1905
            jp = jp + 1
 
1906
 
 
1907
*           **** phase factor does not matter therefore ****
 
1908
*           **** (-i)^l is the same as (i)^l in the     ****
 
1909
*           **** Rayleigh scattering formula            ****
 
1910
 
 
1911
*           *** current function is s or d ****
 
1912
            if (sd_function) then
 
1913
               call Pack_tc_Mul(1,dbl_mb(shift),
 
1914
     >                            dcpl_mb(exi(1)),
 
1915
     >                            dcpl_mb(prjtmp(1)))
 
1916
 
 
1917
*           *** current function is p or f ****
 
1918
            else
 
1919
               call Pack_tc_iMul(1,dbl_mb(shift),
 
1920
     >                            dcpl_mb(exi(1)),
 
1921
     >                            dcpl_mb(prjtmp(1)))
 
1922
 
 
1923
            end if
 
1924
            call Pack_cc_indot(1,nn,
 
1925
     >                       psi1,
 
1926
     >                       dcpl_mb(prjtmp(1)),
 
1927
     >                       dbl_mb(sw1a(1)+(jp-1)*nn))
 
1928
         end do
 
1929
      end do
 
1930
      call D3dB_Vector_SumAll(nn*n_prj_indx,dbl_mb(sw1a(1)))
 
1931
 
 
1932
 
 
1933
*     **** Compute sw2  ****
 
1934
      eh_atom = 0.0d0
 
1935
      do ip=1,nion_prj_indx
 
1936
         ii          = int_mb(ii_prj_indx(1)+ip-1)
 
1937
         ia          = int_mb(ia_prj_indx(1)+ip-1)
 
1938
         nproj       = int_mb(nproj_prj_indx(1)+ip-1)
 
1939
         swstart     = int_mb(swstart_prj_indx(1)+ip-1)
 
1940
 
 
1941
 
 
1942
*        **** sw2 = Gijl*sw1 ******
 
1943
         Gijl_indx = psi_data_get_ptr(int_mb(Gijl(1)+ia-1),1)
 
1944
         call Multiply_Gijl_sw1(nn,
 
1945
     >                          nproj,
 
1946
     >                          int_mb(nmax(1)+ia-1),
 
1947
     >                          int_mb(lmax(1)+ia-1),
 
1948
     >                          int_mb(n_projector(1)
 
1949
     >                                 + (ia-1)*(nmax_max*lmmax_max)),
 
1950
     >                          int_mb(l_projector(1)
 
1951
     >                                 + (ia-1)*(nmax_max*lmmax_max)),
 
1952
     >                          int_mb(m_projector(1)
 
1953
     >                                 + (ia-1)*(nmax_max*lmmax_max)),
 
1954
     >                          dbl_mb(Gijl_indx),
 
1955
     >                          dbl_mb(sw1a(1)+swstart*nn),
 
1956
     >                          dbl_mb(sw2a(1)+swstart*nn))
 
1957
 
 
1958
     
 
1959
*        **** paw operations #1 - generate it's compcharge, add atomic coulomb, and add xc potential ****
 
1960
         if ((int_mb(psp_type(1)+ia-1).eq.4)) then
 
1961
 
 
1962
*           **** paw atom - generate it's atomic density matrix ****
 
1963
            call psp_gen_density_matrix(ispin,ne,nproj,
 
1964
     >                                  dbl_mb(sw1a(1)+swstart*nn),
 
1965
     >                                  dbl_mb(wtmp(1)))
 
1966
 
 
1967
*           **** paw atom - generate it's compcharge ***
 
1968
            call nwpw_compcharge_gen_Qlm(ii,ia,ispin,nproj,
 
1969
     >                                   dbl_mb(wtmp(1)))
 
1970
 
 
1971
*           **** atomic coulomb matrix - sw2 = sw2 + Vhatomijl*sw1  ****
 
1972
            call nwpw_compcharge_coulomb_atom(ii,ia,ispin,ne,nproj,
 
1973
     >                                   dbl_mb(wtmp(1)),
 
1974
     >                                   dbl_mb(sw1a(1)+swstart*nn),
 
1975
     >                                   dbl_mb(sw2a(1)+swstart*nn),
 
1976
     >                                   eh_atom)
 
1977
c            write(*,*) "eh_atom=",ip,ii,eh_atom
 
1978
 
 
1979
*           **** xc matrix - sw2 = sw2 + Vxcijl*sw1 ******
 
1980
            call nwpw_xc_solve(ii,ia,
 
1981
     >        int_mb(n1dgrid(1)+ia-1),
 
1982
     >        int_mb(n1dbasis(1)+ia-1),
 
1983
     >        dbl_mb(psi_data_get_chnk(int_mb(phi_ae(1)+ia-1))),
 
1984
     >        dbl_mb(psi_data_get_chnk(int_mb(phi_ps(1)+ia-1))),
 
1985
     >        dbl_mb(psi_data_get_chnk(int_mb(dphi_ae(1)+ia-1))),
 
1986
     >        dbl_mb(psi_data_get_chnk(int_mb(dphi_ps(1)+ia-1))),
 
1987
     >        dbl_mb(psi_data_get_chnk(int_mb(core_ae(1)+ia-1))),
 
1988
     >        dbl_mb(psi_data_get_chnk(int_mb(core_ps(1)+ia-1))),
 
1989
     >        dbl_mb(psi_data_get_chnk(int_mb(core_ae_prime(1)+ia-1))),
 
1990
     >        dbl_mb(psi_data_get_chnk(int_mb(core_ps_prime(1)+ia-1))),
 
1991
     >        dbl_mb(psi_data_get_chnk(int_mb(rgrid(1)+ia-1))),
 
1992
     >        dbl_mb(log_amesh(1)+ia-1),
 
1993
     >        ispin,ne,nproj,
 
1994
     >        dbl_mb(sw1a(1)+swstart*nn),dbl_mb(sw2a(1)+swstart*nn))
 
1995
 
 
1996
         end if
 
1997
      end do
 
1998
 
 
1999
 
 
2000
*     **** paw operations #2 - generate vcmp,.... and Gaussian Multipole  ****
 
2001
      if (pawexist) then
 
2002
         call D3dB_nx(1,nx)
 
2003
         call D3dB_ny(1,ny)
 
2004
         call D3dB_nz(1,nz)
 
2005
         scal1 = 1.0d0/dble(nx*ny*nz)
 
2006
 
 
2007
         if (periodic) then
 
2008
            call Pack_npack(0,npack0)
 
2009
         else
 
2010
            call D3dB_nfft3d(1,npack0)
 
2011
         end if
 
2012
 
 
2013
*        **** zero out dE/dQlm array ****
 
2014
         call nwpw_compcharge_zero_dE_Qlm()
 
2015
 
 
2016
         if (use_grid_cmp) then
 
2017
            value = MA_push_get(mt_dcpl,npack0,'dng_cmp',
 
2018
     >                          dng_cmp(2),dng_cmp(1))
 
2019
            value = value.and.
 
2020
     >              MA_push_get(mt_dcpl,npack0,'vcmp',
 
2021
     >                          vcmp(2),vcmp(1))
 
2022
            if (.not.value) 
 
2023
     >      call errquit('v_nonlocal:out of stack',4, MA_ERR)
 
2024
 
 
2025
            call nwpw_compcharge_gen_zv_cmp(dbl_mb(zv(1)),
 
2026
     >                                      dcpl_mb(dng_cmp(1)))
 
2027
            call coulomb_v(dcpl_mb(dng_cmp(1)),dcpl_mb(vcmp_tmp(1)))
 
2028
            call Pack_cc_dot(0,dcpl_mb(dng_cmp(1)),
 
2029
     >                         dcpl_mb(vcmp_tmp(1)),ff(1))
 
2030
c            write(*,*) "e_i-i=",0.5d0*ff(1)
 
2031
 
 
2032
            call nwpw_compcharge_gen_dn_cmp(ispin,dcpl_mb(dng_cmp(1)))
 
2033
 
 
2034
c            call nwpw_compcharge_gen_dn_cmp_zv(ispin,dbl_mb(zv(1)),
 
2035
c     >                                         dcpl_mb(dng_cmp(1)))
 
2036
 
 
2037
            call coulomb_v(dcpl_mb(dng_cmp(1)),dcpl_mb(vcmp_tmp(1)))
 
2038
            call Pack_cc_Sum(0,
 
2039
     >                       dcpl_mb(vc_tmp(1)),
 
2040
     >                       dcpl_mb(vcmp_tmp(1)),
 
2041
     >                       dcpl_mb(vcmp(1)))
 
2042
 
 
2043
            pcount = 0
 
2044
            do ip=1,nion_prj_indx
 
2045
               ii          = int_mb(ii_prj_indx(1)+ip-1)
 
2046
               ia          = int_mb(ia_prj_indx(1)+ip-1)
 
2047
               if ((int_mb(psp_type(1)+ia-1).eq.4)) then
 
2048
                  mult_l = nwpw_compcharge_mult_l(ia)
 
2049
                  do l=0,mult_l
 
2050
                  do m=-l,l
 
2051
                     if (mod(pcount,np_j).eq.taskid_j) then
 
2052
                        call nwpw_compcharge_gen_glm(ii,l,m,
 
2053
     >                                      dcpl_mb(dng_cmp(1)))
 
2054
                        call Pack_cc_idot(0,dcpl_mb(dng_cmp(1)),
 
2055
     >                                      dcpl_mb(vcmp(1)),ff(1))
 
2056
c                        write(*,*) "ff=",ff(1)*omega
 
2057
                        call nwpw_compcharge_add_dE_Qlm(ispin,ii,l,m,
 
2058
     >                                                  ff(1)*omega)
 
2059
                     end if
 
2060
                     pcount = pcount + 1
 
2061
                  end do
 
2062
                  end do
 
2063
               end if
 
2064
            end do
 
2065
            value =           MA_pop_stack(vcmp(2))
 
2066
            value = value.and.MA_pop_stack(dng_cmp(2))
 
2067
            if (.not.value) 
 
2068
     >      call errquit('v_nonlocal:popping stack',4, MA_ERR)
 
2069
 
 
2070
         else
 
2071
 
 
2072
            call nwpw_compcharge_add_dEmult_Qlm(ispin)
 
2073
 
 
2074
            value = MA_push_get(mt_dcpl,npack0,'dng_cmp',
 
2075
     >                          dng_cmp(2),dng_cmp(1))
 
2076
            value = value.and.
 
2077
     >              MA_push_get(mt_dcpl,npack0,'dng_cmp_smooth',
 
2078
     >                          dng_cmp_smooth(2),dng_cmp_smooth(1))
 
2079
            value = value.and.
 
2080
     >              MA_push_get(mt_dcpl,npack0,'vcmp',
 
2081
     >                          vcmp(2),vcmp(1))
 
2082
            value = value.and.
 
2083
     >              MA_push_get(mt_dcpl,npack0,'vcmp_smooth',
 
2084
     >                          vcmp_smooth(2),vcmp_smooth(1))
 
2085
            if (.not.value) 
 
2086
     >        call errquit('v_nonlocal:out of stack',4, MA_ERR)
 
2087
 
 
2088
            call nwpw_compcharge_gen_dn_cmp2_zv(ispin,dbl_mb(zv(1)),
 
2089
     >                                       dcpl_mb(dng_cmp(1)),
 
2090
     >                                       dcpl_mb(dng_cmp_smooth(1)))
 
2091
 
 
2092
            !*** compute hartree potential of ntilde + ncmp_tilde ***
 
2093
            !*** compute hartree potential ncmp   - ncmp_tilde ***
 
2094
            if (periodic) then
 
2095
               call coulomb_v(dcpl_mb(dng_cmp(1)),dcpl_mb(vcmp(1)))
 
2096
               call coulomb_v(dcpl_mb(dng_cmp_smooth(1)),
 
2097
     >                        dcpl_mb(vcmp_smooth(1)))
 
2098
 
 
2099
               call Pack_c_Copy(0,dcpl_mb(vcmp(1)),dcpl_mb(vcmp_tmp(1)))
 
2100
               call Pack_cc_Sub(0,dcpl_mb(vcmp_tmp(1)),
 
2101
     >                            dcpl_mb(vcmp_smooth(1)),
 
2102
     >                            dcpl_mb(vcmp(1)))
 
2103
               call Pack_cc_Sum2(0,dcpl_mb(vc_tmp(1)),
 
2104
     >                             dcpl_mb(vcmp_smooth(1)))
 
2105
 
 
2106
 
 
2107
            else
 
2108
c               call Pack_c_unpack(0,dcpl_mb(dng_cmp(1)))
 
2109
c               call D3dB_cr_fft3b(1,dcpl_mb(dng_cmp(1)))
 
2110
c               call coulomb2_v(dcpl_mb(dng_cmp(1)),dcpl_mb(vcmp(1)))
 
2111
c               call D3dB_r_SMul1(1,scal1,dcpl_mb(vcmp(1)))
 
2112
c               call D3dB_rc_fft3f(1,dcpl_mb(vcmp(1)))
 
2113
c               call Pack_c_pack(0,dcpl_mb(vcmp(1)))
 
2114
c            
 
2115
c               call Pack_c_unpack(0,dcpl_mb(dng_cmp_smooth(1)))
 
2116
c               call D3dB_cr_fft3b(1,dcpl_mb(dng_cmp_smooth(1)))
 
2117
c               call coulomb2_v(dcpl_mb(dng_cmp_smooth(1)),
 
2118
c     >                        dcpl_mb(vcmp_smooth(1)))
 
2119
c               call D3dB_r_SMul1(1,scal1,dcpl_mb(vcmp(1)))
 
2120
c               call D3dB_rc_fft3f(1,dcpl_mb(vcmp_smooth(1)))
 
2121
c               call Pack_c_pack(0,dcpl_mb(vcmp_smooth(1)))
 
2122
            end if
 
2123
         
 
2124
 
 
2125
            pcount = 0
 
2126
            do ip=1,nion_prj_indx
 
2127
               ii          = int_mb(ii_prj_indx(1)+ip-1)
 
2128
               ia          = int_mb(ia_prj_indx(1)+ip-1)
 
2129
               if ((int_mb(psp_type(1)+ia-1).eq.4)) then
 
2130
                  mult_l = nwpw_compcharge_mult_l(ia)
 
2131
                  do l=0,mult_l
 
2132
                  do m=-l,l
 
2133
                     if (mod(pcount,np_j).eq.taskid_j) then
 
2134
                        call nwpw_compcharge_gen_glm2(ii,l,m,
 
2135
     >                           dcpl_mb(dng_cmp(1)),
 
2136
     >                           dcpl_mb(dng_cmp_smooth(1)))
 
2137
                        call Pack_cc_idot(0,
 
2138
     >                                    dcpl_mb(dng_cmp(1)),
 
2139
     >                                    dcpl_mb(vcmp_smooth(1)),
 
2140
     >                                    ff(1))
 
2141
                        call Pack_cc_idot(0,
 
2142
     >                                    dcpl_mb(dng_cmp_smooth(1)),
 
2143
     >                                    dcpl_mb(vcmp(1)),
 
2144
     >                                    ff(2))
 
2145
c                        write(*,*) "ff=",ff(1)*omega,ff(2)*omega
 
2146
                        call nwpw_compcharge_add_dE_Qlm(ispin,ii,l,m,
 
2147
     >                                           (ff(1)+ff(2))*omega)
 
2148
                     end if
 
2149
                     pcount = pcount + 1
 
2150
                  end do
 
2151
                  end do
 
2152
               end if
 
2153
            end do
 
2154
 
 
2155
            value =           MA_pop_stack(vcmp_smooth(2))
 
2156
            value = value.and.MA_pop_stack(vcmp(2))
 
2157
            value = value.and.MA_pop_stack(dng_cmp_smooth(2))
 
2158
            value = value.and.MA_pop_stack(dng_cmp(2))
 
2159
            if (.not.value) 
 
2160
     >        call errquit('v_nonlocal:popping stack',5,MA_ERR)
 
2161
         end if
 
2162
 
 
2163
         call nwpw_compcharge_reduce_dE_Qlm()
 
2164
 
 
2165
         do ip=1,nion_prj_indx
 
2166
            ii          = int_mb(ii_prj_indx(1)+ip-1)
 
2167
            ia          = int_mb(ia_prj_indx(1)+ip-1)
 
2168
            nproj       = int_mb(nproj_prj_indx(1)+ip-1)
 
2169
            swstart     = int_mb(swstart_prj_indx(1)+ip-1)
 
2170
            call nwpw_compcharge_gen_sw2(ii,ia,ispin,ne,nproj,
 
2171
     >                                   dbl_mb(sw1a(1)+swstart*nn),
 
2172
     >                                   dbl_mb(sw2a(1)+swstart*nn))
 
2173
         end do
 
2174
     
 
2175
 
 
2176
      end if
 
2177
 
 
2178
 
 
2179
*     **** do Kleinman-Bylander Multiplication ****
 
2180
      call dscal(nn*n_prj_indx,scal,dbl_mb(sw2a(1)),1)
 
2181
 
 
2182
*     **** apply the sw2 to psi ****
 
2183
      jp  = 0
 
2184
      do iip=1,nion_prj_indx,nprj_mult 
 
2185
 
 
2186
         swstart = int_mb(swstart_prj_indx(1)+iip-1)
 
2187
         l1      = 0
 
2188
         iipmax  = (iip+nprj_mult-1)
 
2189
         if (iipmax.gt.nion_prj_indx) iipmax = nion_prj_indx
 
2190
 
 
2191
         do ip=iip,iipmax
 
2192
            ii    = int_mb(ii_prj_indx(1)+ip-1)
 
2193
            ia    = int_mb(ia_prj_indx(1)+ip-1)
 
2194
            nproj = int_mb(nproj_prj_indx(1)+ip-1)
 
2195
 
 
2196
*           **** structure factor and local pseudopotential ****
 
2197
            call strfac_pack(1,ii,dcpl_mb(exi(1)))
 
2198
 
 
2199
            do l=1,nproj
 
2200
               shift       = int_mb(shift_prj_indx(1)+jp)
 
2201
               l_prj       = int_mb(l_prj_prj_indx(1)+jp)
 
2202
               sd_function = log_mb(sd_function_prj_indx(1)+jp)
 
2203
               jp = jp + 1
 
2204
 
 
2205
*              **** phase factor does not matter therefore ****
 
2206
*              **** (-i)^l is the same as (i)^l in the     ****
 
2207
*              **** Rayleigh scattering formula            ****
 
2208
 
 
2209
*              *** current function is s or d ****
 
2210
               if (sd_function) then
 
2211
                  call Pack_tc_Mul(1,dbl_mb(shift),
 
2212
     >                               dcpl_mb(exi(1)),
 
2213
     >                               dcpl_mb(prjtmp(1)+l1*npack1))
 
2214
 
 
2215
*              *** current function is p or f ****
 
2216
               else
 
2217
                  call Pack_tc_iMul(1,dbl_mb(shift),
 
2218
     >                               dcpl_mb(exi(1)),
 
2219
     >                               dcpl_mb(prjtmp(1)+l1*npack1))
 
2220
               end if
 
2221
 
 
2222
*              ***** scale (sw2a) psp by factor - used for generating antiferromagnetic structures ****
 
2223
*              **** nwchem input: pspspin up/down scale l ion_numbers                              ****
 
2224
               if (pspspin) then
 
2225
                  if (log_mb(pspspin_upions(1)+ii-1).and.
 
2226
     >               (l_prj.eq.pspspin_upl))
 
2227
     >            call dscal(ne(1),pspspin_upscale,
 
2228
     >                       dbl_mb(sw2a(1)+(swstart+l-1)*nn),1)
 
2229
                  if (log_mb(pspspin_downions(1)+ii-1).and.
 
2230
     >               (l_prj.eq.pspspin_downl))
 
2231
     >            call dscal(ne(2),pspspin_downscale,
 
2232
     >                       dbl_mb(sw2a(1)+(swstart+l-1)*nn+ne(1)),1)
 
2233
               end if
 
2234
 
 
2235
               l1 = l1 + 1
 
2236
            end do
 
2237
         end do
 
2238
 
 
2239
         call DGEMM('N','T',2*npack1,nn,l1,
 
2240
     >             (-1.0d0),
 
2241
     >             dcpl_mb(prjtmp(1)), 2*npack1,
 
2242
     >             dbl_mb(sw2a(1)+swstart*nn),   nn,
 
2243
     >             (1.0d0),
 
2244
     >             psi2,               2*npack1)
 
2245
 
 
2246
 
 
2247
         if (move) then
 
2248
           l1 = 0
 
2249
           do ip=iip,iipmax
 
2250
              ii    = int_mb(ii_prj_indx(1)+ip-1)
 
2251
              ia    = int_mb(ia_prj_indx(1)+ip-1)
 
2252
              nproj = int_mb(nproj_prj_indx(1)+ip-1)
 
2253
 
 
2254
              do l=1,nproj
 
2255
                 do n=1,nn
 
2256
                    call Pack_cct_iconjgMul(1,
 
2257
     >                                 dcpl_mb(prjtmp(1)+l1*npack1),
 
2258
     >                                 psi1(1+(n-1)*npack1),
 
2259
     >                                 dbl_mb(xtmp(1)))
 
2260
                   call Pack_tt_idot(1,dbl_mb(G(1)),dbl_mb(xtmp(1)),
 
2261
     >                               dbl_mb(sum(1)+3*(n-1)))
 
2262
                   call Pack_tt_idot(1,dbl_mb(G(2)),dbl_mb(xtmp(1)),
 
2263
     >                               dbl_mb(sum(1)+1+3*(n-1)))
 
2264
                   call Pack_tt_idot(1,dbl_mb(G(3)),dbl_mb(xtmp(1)),
 
2265
     >                               dbl_mb(sum(1)+2+3*(n-1)))
 
2266
 
 
2267
                 end do
 
2268
                 call D3dB_Vector_SumAll(3*(nn),dbl_mb(sum(1)))
 
2269
 
 
2270
                 !**** fractional weighting ****
 
2271
                 if (fractional) then
 
2272
                  do n=1,nn
 
2273
                   call Dneall_qton(n,i)
 
2274
                   dbl_mb(sum(1)+3*(n-1))  
 
2275
     >                =dbl_mb(sum(1)  +3*(n-1))*occ(i)
 
2276
                   dbl_mb(sum(1)+1+3*(n-1))
 
2277
     >                =dbl_mb(sum(1)+1+3*(n-1))*occ(i)
 
2278
                   dbl_mb(sum(1)+2+3*(n-1))
 
2279
     >                =dbl_mb(sum(1)+2+3*(n-1))*occ(i)
 
2280
                  end do
 
2281
                 end if
 
2282
 
 
2283
c                 ff(1) = 0.0d0
 
2284
c                 ff(2) = 0.0d0
 
2285
c                 ff(3) = 0.0d0
 
2286
c                 do n=1,nn
 
2287
c                    ff(1) = ff(1) + 2.0d0*dbl_mb(sw2a(1)+swstart*nn+n-1+(l-1)*nn) !// change
 
2288
c     >                                   *dbl_mb(sum(1)+  3*(n-1))
 
2289
c                    ff(2) = ff(2) + 2.0d0*dbl_mb(sw2a(1)+swstart*nn+n-1+(l-1)*nn) !// change
 
2290
c     >                                   *dbl_mb(sum(1)+1+3*(n-1))
 
2291
c                    ff(3) = ff(3) + 2.0d0*dbl_mb(sw2a(1)+swstart*nn+n-1+(l-1)*nn) !// change
 
2292
c     >                                   *dbl_mb(sum(1)+2+3*(n-1))
 
2293
c                 end do
 
2294
                 ff(1) =2.0d0*ddot(nn,dbl_mb(sw2a(1)+(swstart+l1)*nn),1,
 
2295
     >                                dbl_mb(sum(1)),3)
 
2296
                 ff(2) =2.0d0*ddot(nn,dbl_mb(sw2a(1)+(swstart+l1)*nn),1,
 
2297
     >                                dbl_mb(sum(1)+1),3)
 
2298
                 ff(3) =2.0d0*ddot(nn,dbl_mb(sw2a(1)+(swstart+l1)*nn),1,
 
2299
     >                                dbl_mb(sum(1)+2),3)
 
2300
                 call D1dB_Vector_SumAll(3,ff)
 
2301
                 fion(1,ii) = fion(1,ii)  + ff(1)*(3-ispin)
 
2302
                 fion(2,ii) = fion(2,ii)  + ff(2)*(3-ispin)
 
2303
                 fion(3,ii) = fion(3,ii)  + ff(3)*(3-ispin)
 
2304
 
 
2305
                 l1 = l1 + 1
 
2306
              end do !** l **
 
2307
           end do !** ip **
 
2308
         end if !** move **
 
2309
 
 
2310
      end do
 
2311
 
 
2312
      value = .true.
 
2313
      if (move) then
 
2314
      value = value.and.MA_pop_stack(sum(2))
 
2315
      value = value.and.MA_pop_stack(xtmp(2))
 
2316
      end if
 
2317
      value = value.and.MA_pop_stack(exi(2))
 
2318
      if (.not.value) call errquit('v_nonlocal: popping stack',3,
 
2319
     &       MA_ERR)
 
2320
 
 
2321
      call nwpw_timing_end(6)
 
2322
      return 
 
2323
      end
 
2324
 
 
2325
 
 
2326
 
 
2327
 
1509
2328
 
1510
2329
 
1511
2330
      subroutine Multiply_Gijl_sw1(nn,nprj,nmax,lmax,
1546
2365
      return
1547
2366
      end
1548
2367
 
1549
 
 
1550
 
*     ***********************************
1551
 
*     *                                 *
1552
 
*     *            f_vnonlocal          *
1553
 
*     *                                 *
1554
 
*     ***********************************
1555
 
 
1556
 
      subroutine f_vnonlocal(ispin,ne,psi1,fion,fractional,occ)
 
2368
*     ***********************************
 
2369
*     *                                 *
 
2370
*     *            E_vnonlocal          *
 
2371
*     *                                 *
 
2372
*     ***********************************
 
2373
      real*8 function E_vnonlocal(ispin,ne,fractional,occ)
 
2374
      implicit none
 
2375
      integer ispin,ne(2)
 
2376
      logical fractional
 
2377
      real*8 occ(*)
 
2378
 
 
2379
#include "mafdecls.fh"
 
2380
#include "psp.fh"
 
2381
#include "errquit.fh"
 
2382
 
 
2383
*     **** local variables ****
 
2384
      real*8 E
 
2385
 
 
2386
*     **** external functions ****
 
2387
      real*8   E_vnonlocal_sub
 
2388
      external E_vnonlocal_sub
 
2389
 
 
2390
      E = E_vnonlocal_sub(ne(1)+ne(2),n_prj_indx,
 
2391
     >                    dbl_mb(sw1a(1)),dbl_mb(sw2a(1)),
 
2392
     >                    fractional,occ)
 
2393
      call D1dB_SumAll(E)
 
2394
      if (ispin.eq.1) E = E+E
 
2395
 
 
2396
      E_vnonlocal = E
 
2397
      return
 
2398
      end
 
2399
 
 
2400
      real*8 function E_vnonlocal_sub(nn,nprjs,sw1,sw2,fractional,occ)
 
2401
      implicit none
 
2402
      integer nn,nprjs
 
2403
      real*8 sw1(nn,nprjs),sw2(nn,nprjs)
 
2404
      logical fractional
 
2405
      real*8 occ(*)
 
2406
 
 
2407
*     **** local variables ****
 
2408
      integer i,ip,n
 
2409
      real*8 sum
 
2410
 
 
2411
      sum = 0.0d0
 
2412
      if (fractional) then
 
2413
         do ip=1,nprjs
 
2414
            do n=1,nn
 
2415
               call Dneall_qton(n,i)
 
2416
               sum = sum + sw1(n,ip)*sw2(n,ip)*occ(i)
 
2417
            end do
 
2418
         end do
 
2419
      else
 
2420
         do ip=1,nprjs
 
2421
            do n=1,nn
 
2422
               sum = sum + sw1(n,ip)*sw2(n,ip)
 
2423
            end do
 
2424
         end do
 
2425
      end if
 
2426
 
 
2427
      E_vnonlocal_sub = sum
 
2428
      return
 
2429
      end
 
2430
 
 
2431
*     ***********************************
 
2432
*     *                                 *
 
2433
*     *            f_vnonlocal_old      *
 
2434
*     *                                 *
 
2435
*     ***********************************
 
2436
 
 
2437
      subroutine f_vnonlocal_old(ispin,ne,psi1,fion,fractional,occ)
1557
2438
      implicit none
1558
2439
      integer    ispin,ne(2)
1559
2440
      complex*16 psi1(*)
1803
2684
      end
1804
2685
 
1805
2686
 
 
2687
 
 
2688
*     ***********************************
 
2689
*     *                                 *
 
2690
*     *        f_vnonlocal              *
 
2691
*     *                                 *
 
2692
*     ***********************************
 
2693
 
 
2694
*    This routine computes the Kleinman-Bylander non-local 
 
2695
* pseudopotential projection.
 
2696
*
 
2697
*  Note - This routine was restructured 12-1-2013 to handle PAW operators.
 
2698
*
 
2699
*  To Do -  For very large numbers of atoms the code will need to be restructured
 
2700
*           to distribute the sw1a and sw2a matrices over np_i.  Basically, if the orthogonalization matrices
 
2701
*           need to be distributed then sw1a and sw2a will need to be distributed as well.  A simple algorithm
 
2702
*           to do this will be to keep the loop structure the same but instead of the
 
2703
*           call to D3dB_Vector_SumAll(nn*n_prj_indx,dbl_mb(sw1a(1))), this will need to be, 
 
2704
*           changed to D3dB_Vector_SumAll(nn*nproj,dbl_mb(sw1(1))) and then place sw1 --> sw1a(distributed)
 
2705
*
 
2706
      subroutine f_vnonlocal(ispin,ne,psi1,fion,fractional,occ)
 
2707
      implicit none
 
2708
      integer    ispin,ne(2)
 
2709
      complex*16 psi1(*)
 
2710
      real*8 fion(3,*)
 
2711
      logical fractional
 
2712
      real*8 occ(*)
 
2713
 
 
2714
#include "mafdecls.fh"
 
2715
#include "psp.fh"
 
2716
#include "errquit.fh"
 
2717
 
 
2718
*     *** local variables ***
 
2719
      integer G(3),npack1,nion,nu
 
2720
      integer i,j,ii,ia,l,n,nn,l1,ip,iip,iipmax,jp,swstart
 
2721
      integer k,shift,l_prj,nproj,Gijl_indx
 
2722
      real*8  omega,scal,ff(3)
 
2723
      complex*16 ctmp
 
2724
      integer exi(2),xtmp(2),sum(2)
 
2725
      logical value,sd_function
 
2726
      real*8 vmm(50)
 
2727
      integer ld_ptr
 
2728
 
 
2729
*     **** external functions ****
 
2730
      logical  is_sORd
 
2731
      integer  ion_nion,ion_katm,Pack_G_indx
 
2732
      integer  psi_data_get_ptr,psi_data_get_chnk
 
2733
      real*8   lattice_omega,ddot
 
2734
      external is_sORd
 
2735
      external ion_nion,ion_katm,Pack_G_indx
 
2736
      external psi_data_get_ptr,psi_data_get_chnk
 
2737
      external lattice_omega,ddot
 
2738
 
 
2739
      call nwpw_timing_start(6) 
 
2740
 
 
2741
 
 
2742
*     **** allocate local memory ****
 
2743
      nion = ion_nion()
 
2744
      nn = ne(1)+ne(2)
 
2745
      call Pack_npack(1,npack1)
 
2746
 
 
2747
      call psp_prj_indx_alloc_sw1a_sw2a(nn)
 
2748
      value = MA_push_get(mt_dcpl,npack1,'exi',exi(2), exi(1))
 
2749
      value = value.and.MA_push_get(mt_dbl,npack1,
 
2750
     >                               'xtmp',xtmp(2),xtmp(1))
 
2751
      value = value.and.MA_push_get(mt_dbl,3*nn,'sum',sum(2),sum(1))
 
2752
      if (.not.value) 
 
2753
     >  call errquit('f_vnonlocal:out of stack',0, MA_ERR)
 
2754
 
 
2755
      G(1)  = Pack_G_indx(1,1)
 
2756
      G(2)  = Pack_G_indx(1,2)
 
2757
      G(3)  = Pack_G_indx(1,3)
 
2758
 
 
2759
      omega = lattice_omega()
 
2760
      scal = 1.0d0/(omega)
 
2761
 
 
2762
      jp = 0
 
2763
      do ip=1,nion_prj_indx
 
2764
         ii          = int_mb(ii_prj_indx(1)+ip-1)
 
2765
         ia          = int_mb(ia_prj_indx(1)+ip-1)
 
2766
         nproj       = int_mb(nproj_prj_indx(1)+ip-1)
 
2767
 
 
2768
*        **** structure factor and local pseudopotential ****
 
2769
         call strfac_pack(1,ii,dcpl_mb(exi(1)))
 
2770
 
 
2771
         do l=1,nproj
 
2772
            shift       = int_mb(shift_prj_indx(1)+jp)
 
2773
            sd_function = log_mb(sd_function_prj_indx(1)+jp)
 
2774
            jp = jp + 1
 
2775
 
 
2776
*           **** phase factor does not matter therefore ****
 
2777
*           **** (-i)^l is the same as (i)^l in the     ****
 
2778
*           **** Rayleigh scattering formula            ****
 
2779
 
 
2780
*           *** current function is s or d ****
 
2781
            if (sd_function) then
 
2782
               call Pack_tc_Mul(1,dbl_mb(shift),
 
2783
     >                            dcpl_mb(exi(1)),
 
2784
     >                            dcpl_mb(prjtmp(1)))
 
2785
 
 
2786
*           *** current function is p or f ****
 
2787
            else
 
2788
               call Pack_tc_iMul(1,dbl_mb(shift),
 
2789
     >                            dcpl_mb(exi(1)),
 
2790
     >                            dcpl_mb(prjtmp(1)))
 
2791
 
 
2792
            end if
 
2793
            call Pack_cc_indot(1,nn,
 
2794
     >                       psi1,
 
2795
     >                       dcpl_mb(prjtmp(1)),
 
2796
     >                       dbl_mb(sw1a(1)+(jp-1)*nn))
 
2797
         end do
 
2798
      end do
 
2799
      call D3dB_Vector_SumAll(nn*n_prj_indx,dbl_mb(sw1a(1)))
 
2800
 
 
2801
 
 
2802
*     **** Compute sw2  ****
 
2803
      do ip=1,nion_prj_indx
 
2804
         ii          = int_mb(ii_prj_indx(1)+ip-1)
 
2805
         ia          = int_mb(ia_prj_indx(1)+ip-1)
 
2806
         nproj       = int_mb(nproj_prj_indx(1)+ip-1)
 
2807
         swstart     = int_mb(swstart_prj_indx(1)+ip-1)
 
2808
 
 
2809
*        **** sw2 = Gijl*sw1 ******
 
2810
         Gijl_indx = psi_data_get_ptr(int_mb(Gijl(1)+ia-1),1)
 
2811
         call Multiply_Gijl_sw1(nn,
 
2812
     >                          nproj,
 
2813
     >                          int_mb(nmax(1)+ia-1),
 
2814
     >                          int_mb(lmax(1)+ia-1),
 
2815
     >                          int_mb(n_projector(1)
 
2816
     >                                 + (ia-1)*(nmax_max*lmmax_max)),
 
2817
     >                          int_mb(l_projector(1)
 
2818
     >                                 + (ia-1)*(nmax_max*lmmax_max)),
 
2819
     >                          int_mb(m_projector(1)
 
2820
     >                                 + (ia-1)*(nmax_max*lmmax_max)),
 
2821
     >                          dbl_mb(Gijl_indx),
 
2822
     >                          dbl_mb(sw1a(1)+swstart*nn),
 
2823
     >                          dbl_mb(sw2a(1)+swstart*nn))
 
2824
 
 
2825
      end do
 
2826
 
 
2827
*     **** do Kleinman-Bylander Multiplication ****
 
2828
      call dscal(nn*n_prj_indx,scal,dbl_mb(sw2a(1)),1)
 
2829
 
 
2830
*     **** apply the sw2 to psi ****
 
2831
      jp  = 0
 
2832
      do iip=1,nion_prj_indx,nprj_mult 
 
2833
 
 
2834
         swstart = int_mb(swstart_prj_indx(1)+iip-1)
 
2835
         l1      = 0
 
2836
         iipmax  = (iip+nprj_mult-1)
 
2837
         if (iipmax.gt.nion_prj_indx) iipmax = nion_prj_indx
 
2838
 
 
2839
         do ip=iip,iipmax
 
2840
            ii    = int_mb(ii_prj_indx(1)+ip-1)
 
2841
            ia    = int_mb(ia_prj_indx(1)+ip-1)
 
2842
            nproj = int_mb(nproj_prj_indx(1)+ip-1)
 
2843
 
 
2844
*           **** structure factor and local pseudopotential ****
 
2845
            call strfac_pack(1,ii,dcpl_mb(exi(1)))
 
2846
 
 
2847
            do l=1,nproj
 
2848
               shift       = int_mb(shift_prj_indx(1)+jp)
 
2849
               l_prj       = int_mb(l_prj_prj_indx(1)+jp)
 
2850
               sd_function = log_mb(sd_function_prj_indx(1)+jp)
 
2851
               jp = jp + 1
 
2852
 
 
2853
*              **** phase factor does not matter therefore ****
 
2854
*              **** (-i)^l is the same as (i)^l in the     ****
 
2855
*              **** Rayleigh scattering formula            ****
 
2856
 
 
2857
*              *** current function is s or d ****
 
2858
               if (sd_function) then
 
2859
                  call Pack_tc_Mul(1,dbl_mb(shift),
 
2860
     >                               dcpl_mb(exi(1)),
 
2861
     >                               dcpl_mb(prjtmp(1)+l1*npack1))
 
2862
 
 
2863
*              *** current function is p or f ****
 
2864
               else
 
2865
                  call Pack_tc_iMul(1,dbl_mb(shift),
 
2866
     >                               dcpl_mb(exi(1)),
 
2867
     >                               dcpl_mb(prjtmp(1)+l1*npack1))
 
2868
               end if
 
2869
 
 
2870
*              ***** scale (sw2a) psp by factor - used for generating antiferromagnetic structures ****
 
2871
*              **** nwchem input: pspspin up/down scale l ion_numbers                              ****
 
2872
               if (pspspin) then
 
2873
                  if (log_mb(pspspin_upions(1)+ii-1).and.
 
2874
     >               (l_prj.eq.pspspin_upl))
 
2875
     >            call dscal(ne(1),pspspin_upscale,
 
2876
     >                       dbl_mb(sw2a(1)+(swstart+l-1)*nn),1)
 
2877
                  if (log_mb(pspspin_downions(1)+ii-1).and.
 
2878
     >               (l_prj.eq.pspspin_downl))
 
2879
     >            call dscal(ne(2),pspspin_downscale,
 
2880
     >                       dbl_mb(sw2a(1)+(swstart+l-1)*nn+ne(1)),1)
 
2881
               end if
 
2882
 
 
2883
               l1 = l1 + 1
 
2884
            end do
 
2885
         end do
 
2886
 
 
2887
         l1 = 0
 
2888
         do ip=iip,iipmax
 
2889
            ii    = int_mb(ii_prj_indx(1)+ip-1)
 
2890
            ia    = int_mb(ia_prj_indx(1)+ip-1)
 
2891
            nproj = int_mb(nproj_prj_indx(1)+ip-1)
 
2892
 
 
2893
            do l=1,nproj
 
2894
               do n=1,nn
 
2895
                  call Pack_cct_iconjgMul(1,
 
2896
     >                               dcpl_mb(prjtmp(1)+l1*npack1),
 
2897
     >                               psi1(1+(n-1)*npack1),
 
2898
     >                               dbl_mb(xtmp(1)))
 
2899
                 call Pack_tt_idot(1,dbl_mb(G(1)),dbl_mb(xtmp(1)),
 
2900
     >                             dbl_mb(sum(1)+3*(n-1)))
 
2901
                 call Pack_tt_idot(1,dbl_mb(G(2)),dbl_mb(xtmp(1)),
 
2902
     >                             dbl_mb(sum(1)+1+3*(n-1)))
 
2903
                 call Pack_tt_idot(1,dbl_mb(G(3)),dbl_mb(xtmp(1)),
 
2904
     >                             dbl_mb(sum(1)+2+3*(n-1)))
 
2905
 
 
2906
               end do
 
2907
               call D3dB_Vector_SumAll(3*(nn),dbl_mb(sum(1)))
 
2908
 
 
2909
               !**** fractional weighting ****
 
2910
               if (fractional) then
 
2911
                do n=1,nn
 
2912
                 call Dneall_qton(n,i)
 
2913
                 dbl_mb(sum(1)+3*(n-1))  
 
2914
     >              =dbl_mb(sum(1)  +3*(n-1))*occ(i)
 
2915
                 dbl_mb(sum(1)+1+3*(n-1))
 
2916
     >              =dbl_mb(sum(1)+1+3*(n-1))*occ(i)
 
2917
                 dbl_mb(sum(1)+2+3*(n-1))
 
2918
     >              =dbl_mb(sum(1)+2+3*(n-1))*occ(i)
 
2919
                end do
 
2920
               end if
 
2921
 
 
2922
               ff(1) =2.0d0*ddot(nn,dbl_mb(sw2a(1)+(swstart+l1)*nn),1,
 
2923
     >                              dbl_mb(sum(1)),3)
 
2924
               ff(2) =2.0d0*ddot(nn,dbl_mb(sw2a(1)+(swstart+l1)*nn),1,
 
2925
     >                              dbl_mb(sum(1)+1),3)
 
2926
               ff(3) =2.0d0*ddot(nn,dbl_mb(sw2a(1)+(swstart+l1)*nn),1,
 
2927
     >                              dbl_mb(sum(1)+2),3)
 
2928
               call D1dB_Vector_SumAll(3,ff)
 
2929
               fion(1,ii) = fion(1,ii)  + ff(1)*(3-ispin)
 
2930
               fion(2,ii) = fion(2,ii)  + ff(2)*(3-ispin)
 
2931
               fion(3,ii) = fion(3,ii)  + ff(3)*(3-ispin)
 
2932
 
 
2933
               l1 = l1 + 1
 
2934
            end do !** l **
 
2935
         end do !** ip **
 
2936
 
 
2937
      end do !** iip **
 
2938
 
 
2939
      value =           MA_pop_stack(sum(2))
 
2940
      value = value.and.MA_pop_stack(xtmp(2))
 
2941
      value = value.and.MA_pop_stack(exi(2))
 
2942
      if (.not.value) call errquit('f_vnonlocal: popping stack',3,
 
2943
     &       MA_ERR)
 
2944
 
 
2945
      call nwpw_timing_end(6)
 
2946
      return 
 
2947
      end
 
2948
 
 
2949
 
 
2950
 
1806
2951
*     ********************************************
1807
2952
*     *                                          *
1808
2953
*     *             psp_read                     *
1898
3043
      call Parallel_taskid(taskid)
1899
3044
      call Parallel2d_taskid_i(taskid_i)
1900
3045
      rlocal = 0.0d0
 
3046
      n1dgrid = 0
 
3047
      n1dbasis = 0
 
3048
      sigma = 0.0d0
1901
3049
 
1902
3050
      pio = control_parallel_io()
1903
3051
      if (pio) then
1974
3122
         n1dbasis = lmax+1
1975
3123
      end if
1976
3124
 
 
3125
 
1977
3126
      msglen=lmax+1
1978
3127
      call Parallela_Brdcst_values(com_p,MASTER,msglen,rc)
1979
3128
      msglen = 1
2023
3172
*     ***** Miscellaneous paw energies and 1d wavefunctions****
2024
3173
      if (psp_type.eq.4) then
2025
3174
 
2026
 
         !**** allocate and reading hartree_matrix(n1dbasis,n1dbasis,n1dbasis,n1dbasi,2*lmax+1) ****
 
3175
         !**** allocate and reading hartree_matrix(n1dbasis,n1dbasis,n1dbasis,n1dbasis,2*lmax+1) ****
2027
3176
         msglen=n1dbasis*n1dbasis*n1dbasis*n1dbasis
2028
3177
         hartree_tag  = psi_data_alloc(2*lmax+1,msglen)
2029
3178
         indx = psi_data_get_chnk(hartree_tag)
2043
3192
         call Parallela_Brdcst_values(com_p,MASTER,(2*lmax+1)*msglen,
2044
3193
     >                                dbl_mb(indx))
2045
3194
 
2046
 
         !**** allocate and reading comp_pot_matrix(n1dbasis,n1dbasi,2*lmax+1) ****
 
3195
         !**** allocate and reading comp_pot_matrix(n1dbasis,n1dbasis,2*lmax+1) ****
2047
3196
         msglen=n1dbasis*n1dbasis
2048
3197
         comp_pot_tag  = psi_data_alloc(2*lmax+1,msglen)
2049
3198
         indx = psi_data_get_chnk(comp_pot_tag)
2080
3229
        call Parallela_Brdcst_values(com_p,MASTER,msglen,sigma)
2081
3230
        call Parallela_Brdcst_values(com_p,MASTER,msglen,zion)
2082
3231
 
 
3232
 
2083
3233
        eig_tag     = psi_data_alloc(n1dbasis,1)
2084
3234
        phi_ae_tag  = psi_data_alloc(n1dbasis,n1dgrid)
2085
3235
        dphi_ae_tag = psi_data_alloc(n1dbasis,n1dgrid)
2216
3366
 
2217
3367
         !**** allocate and reading r3_matrix(n1dbasis,n1dbasis) ****
2218
3368
         msglen=n1dbasis*n1dbasis
2219
 
c         write(*,*) "n1dbasis,msglen=",n1dbasis,msglen
2220
3369
         r3_tag = psi_data_alloc(1,msglen)
2221
3370
         indx   = psi_data_get_chnk(r3_tag)
2222
3371
         if (taskid_p.eq.MASTER) then
2291
3440
      external     control_boundry
2292
3441
      external     ion_atom
2293
3442
      
 
3443
      call nwpw_timing_start(50)
2294
3444
 
2295
3445
      call D3dB_nfft3d(1,nfft3d)
2296
3446
      call Pack_npack(1,npack1)
2461
3611
     > call errquit('psp_readall:error popping stack',0,MA_ERR)
2462
3612
 
2463
3613
 
 
3614
 
2464
3615
*     **** done reading set nprj_max and prjtmp for nonlocal psp operator ****
2465
3616
      call psp_proj_init()
2466
3617
 
 
3618
*     **** initialize prj_indx, used by v_nonlocal_new ****
 
3619
      call psp_prj_indx_init()
 
3620
 
 
3621
 
2467
3622
      if (pawexist) then
2468
3623
         call nwpw_compcharge_init(ion_nion(),npsp,
2469
3624
     >                             int_mb(nprj(1)),
2475
3630
     >                             int_mb(l_projector(1)),
2476
3631
     >                             int_mb(m_projector(1)),
2477
3632
     >                             int_mb(b_projector(1)),
2478
 
     >                             int_mb(comp_charge_matrix(1)))
 
3633
     >                             int_mb(comp_charge_matrix(1)),
 
3634
     >                             int_mb(hartree_matrix(1)))
2479
3635
 
2480
3636
         call nwpw_xc_init(ion_nion(),npsp,
2481
3637
     >                     int_mb(nprj(1)),
2491
3647
      end if
2492
3648
 
2493
3649
 
 
3650
      call nwpw_timing_end(50)
2494
3651
      return
2495
3652
      end
2496
3653