← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/gradients/grad1.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/gradients/grad1.F
2
2
     $                   frc_kin, frc_wgh, g_force,
3
3
     $                   g_dens, g_wdens, basis, geom, nproc, nat, 
4
4
     $                   max_at_bf, rtdb, oskel )
5
 
c$Id: grad1.F 19708 2010-10-29 18:04:21Z d3y133 $
 
5
c$Id: grad1.F 23020 2012-10-30 01:03:15Z d3y133 $
6
6
 
7
7
C     one electron contribution to RHF, ROHF and UHF gradients
8
8
C     now also UMP2
141
141
 
142
142
C               1el. derivatives
143
143
                if(.not.dobq) then
144
 
                call intd_1eh1(basis,ish1,basis,ish2,lscr,scr,
145
 
     &               lbuf,H)
 
144
                  call intd_1eh1(basis,ish1,basis,ish2,lscr,scr,
 
145
     &                 lbuf,H)
146
146
                else
147
 
                call intd_1epot(basis,ish1,basis,ish2,lscr,scr,
148
 
     &               lbuf,H)
 
147
                  call intd_1epot(basis,ish1,basis,ish2,lscr,scr,
 
148
     &                 lbuf,H)
149
149
                end if
150
 
                  
151
150
 
152
151
C     D x H
153
152
 
250
249
 
251
250
      return
252
251
      end
 
252
c
 
253
C> \brief calculate the gradient terms due to the interaction with the 
 
254
C> COSMO charges
253
255
      subroutine grad_hnd_cos ( H, lbuf, scr, lscr, 
254
 
     $                   dens, wdens, frc_nuc,
255
 
     $                   frc_kin, frc_wgh, g_force,
256
 
     $                   g_dens, g_wdens, basis, geom, nproc, nat, 
 
256
     $                   dens, frc_cos_nucq, frc_cos_elq,
 
257
     $                   frc_cos_qq, 
 
258
     $                   g_dens, basis, geom, nproc, nat, 
257
259
     $                   max_at_bf, rtdb, oskel )
258
 
c$Id: grad1.F 19708 2010-10-29 18:04:21Z d3y133 $
 
260
c$Id: grad1.F 23020 2012-10-30 01:03:15Z d3y133 $
259
261
 
260
 
C     one electron contribution to RHF, ROHF and UHF gradients
261
 
C     now also UMP2
 
262
C     COSMO one electron contribution to RHF, ROHF and UHF gradients
 
263
C     now also UMP2 ??? unlikely as that requires solutions to the
 
264
c     CPHF equation???
 
265
c
 
266
c     Terms included in this subroutine are:
 
267
c     1. Electron - COSMO charge interactions
 
268
c     2. Nuclear - COSMO charge interactions
 
269
c     3. COSMO charge - COSMO charge interactions
 
270
c
 
271
c     Terms NOT included are:
 
272
c     1. All regular QM derivatives
262
273
 
263
274
      implicit none
264
275
 
 
276
#include "nwc_const.fh"
265
277
#include "mafdecls.fh"
 
278
#include "errquit.fh"
266
279
#include "global.fh"
 
280
#include "geomP.fh"
267
281
#include "geom.fh"
 
282
#include "bq.fh"
268
283
#include "bas.fh"
269
284
#include "rtdb.fh"
270
285
#include "sym.fh"
271
286
#include "stdio.fh"
 
287
#include "prop.fh"
272
288
 
273
289
C-------------------------parameters--------------------------------
274
 
      integer lbuf, lscr,
275
 
     $     g_dens,        ! density matrix (summed if ROHF, UHF)
276
 
     $     g_wdens,       ! weighted density (Lagrangian)
277
 
     $     g_force,       ! global force array
278
 
     $     basis, geom, nproc, nat, max_at_bf, rtdb
279
 
 
280
 
      double precision H, ! integral derivatives
281
 
     $     scr, 
282
 
     $     dens,          ! local density block
283
 
     $     wdens,         ! local weighted density block
284
 
     $     frc_nuc, frc_kin, frc_wgh   ! forces arrays
285
 
 
286
 
      dimension H ( lbuf ), frc_nuc(3, nat), frc_kin(3, nat),
287
 
     $          frc_wgh(3, nat), scr(lscr),
288
 
     $          dens(max_at_bf,max_at_bf), wdens(max_at_bf,max_at_bf)
 
290
      integer g_dens !< [Input] the total electron density matrix GA
 
291
      integer basis  !< [Input] the basis set handle
 
292
      integer geom   !< [Input] the geometry handle
 
293
      integer rtdb   !< [Input] the RTDB handle
 
294
      integer lbuf   !< [Input] the length of the integral buffer
 
295
      integer lscr,  !< [Input] the length of the scratch space
 
296
     $     nproc, nat, max_at_bf
 
297
 
 
298
      double precision frc_cos_nucq(3,nat) !< [Output] the forces due
 
299
                                           !< nuclear-COSMO charge 
 
300
                                           !< interaction
 
301
      double precision frc_cos_elq(3,nat)  !< [Output] the forces due
 
302
                                           !< electron-COSMO charge 
 
303
                                           !< interaction
 
304
      double precision frc_cos_qq(3,nat)   !< [Output] the forces due
 
305
                                           !< COSMO charge-COSMO charge 
 
306
                                           !< interaction
 
307
      double precision H(lbuf)   !< [Scratch] the derivative integrals
 
308
      double precision scr(lscr) !< [Scratch] scratch space
 
309
      double precision dens(max_at_bf,max_at_bf) !< [Scratch] local
 
310
                                                 !< density block
289
311
 
290
312
      logical oskel   ! symmetry?
 
313
c
 
314
      double precision dielec, screen, rsolv
 
315
      common/hnd_cospar/dielec, screen ,rsolv
291
316
 
292
317
C-------------------------local variables--------------------------
293
318
 
295
320
     $     iab1f, iab1l, iab2f, iab2l, iac1f, iac1l, iac2f, iac2l,
296
321
     $     if1, il1, if2, il2,
297
322
     $     icart, ic, nint, ip1, ip2
298
 
 
299
 
      double precision dE, qfac
 
323
      integer im1, im2, nprim, ngen, sphcart, ityp1, ityp2
 
324
      integer ich1, ich2
 
325
 
 
326
      integer nefc        ! the number of COSMO charges
 
327
      integer l_efciat    ! the handle of the COSMO charge-atom map
 
328
      integer k_efciat    ! the index of the COSMO charge-atom map
 
329
      integer nefcl       ! the number of COSMO charge for a given atom
 
330
      integer iefc        ! counter over COSMO charges
 
331
      integer iefc_c      ! memory index for COSMO charge coordinates
 
332
      integer iefc_q      ! memory index for COSMO charges
 
333
 
 
334
      double precision dE, qfac, fact, dx, dy, dz, rr
 
335
      double precision invscreen
300
336
 
301
337
      logical status, pointforce
302
338
 
306
342
c     ---- -cosmo- gradient term -----
307
343
      logical odbug
308
344
 
309
 
      odbug=.true.
 
345
      odbug=.false.
310
346
      if(odbug) then
311
347
         write(Luout,*) 'in -grad1_hnd_cos- ...'
312
348
      endif
314
350
 
315
351
      task_size = 1
316
352
      status = rtdb_parallel(.true.) ! Broadcast reads to all processes
317
 
 
 
353
c
318
354
      pointforce = geom_include_bqbq(geom)
319
 
 
 
355
      if (.not.bq_create('cosmo efc bq',cosmo_bq_efc))
 
356
     $   call errquit("grad_hnd_cos: bq_create on cosmo failed",
 
357
     $                0,GEOM_ERR)
 
358
      if (.not.bq_rtdb_load(rtdb,cosmo_bq_efc))
 
359
     $   call errquit('grad_hnd_cos: rtdb load failed for Bq',916,
 
360
     $                rtdb_err)
 
361
      if (.not.bq_ncenter(cosmo_bq_efc,nefc))
 
362
     $   call errquit('grad_hnd_cos: could not retrieve nefc',917,
 
363
     $                GEOM_ERR)
 
364
      if (.not.bq_index_coord(cosmo_bq_efc,iefc_c))
 
365
     $   call errquit('grad_hnd_cos: could not get coordinate index Bq',
 
366
     $                cosmo_bq_efc,MA_ERR)
 
367
      if (.not.bq_index_charge(cosmo_bq_efc,iefc_q))
 
368
     $   call errquit('grad_hnd_cos: could not get charge index Bq',
 
369
     $                cosmo_bq_efc,MA_ERR)
 
370
c
 
371
      if (.not.ma_push_get(MT_INT,nefc,"efciat",l_efciat,k_efciat))
 
372
     $  call errquit("grad_hnd_cos: could not allocate efciat",
 
373
     $               ma_sizeof(MT_BYTE,nefc,MT_INT),MA_ERR)
 
374
      if(.not.rtdb_get(rtdb,'cosmo:efciat',mt_int,nefc,
 
375
     $                 int_mb(k_efciat)))
 
376
     $   call errquit('grad_hnd_cos: rtdb get failed for iatefc',915,
 
377
     $                rtdb_err)
 
378
c
320
379
      call hf_print_set(1)
321
380
 
322
381
      ijatom = -1
343
402
            else
344
403
               qfac = 1.0d0
345
404
            endif
 
405
            qfac = -qfac
346
406
 
347
407
            status = bas_ce2cnr(basis,iat1,iac1f,iac1l)
348
408
            status = bas_ce2cnr(basis,iat2,iac2f,iac2l)
349
409
 
350
410
            call ga_get (g_dens, iab1f,iab1l,iab2f,iab2l,dens,max_at_bf)
351
 
            call ga_get(g_wdens,iab1f,iab1l,iab2f,iab2l,wdens,max_at_bf)
352
411
 
353
412
            do 70, ish1 = iac1f, iac1l
354
413
              if ( iat1.eq.iat2 ) iac2l = ish1
358
417
                status = bas_cn2bfr(basis,ish1,if1,il1)
359
418
                if1 = if1 - iab1f + 1
360
419
                il1 = il1 - iab1f + 1
 
420
c
 
421
c               Work out the number of Cartesian basis functions
 
422
c               The integrals are evaluated in the Cartesian basis set
 
423
c               and then transformed to spherical harmonics. So the
 
424
c               buffer size depends on the number of Cartesian functions
 
425
c
 
426
                status = bas_continfo(basis,ish1,ityp1,nprim,ngen,
 
427
     +                                sphcart)
 
428
                if (sphcart.eq.1.and.ityp1.ge.2) then
 
429
                  im1 = if1 + (ityp1+1)*(ityp1+2)/2 - 1
 
430
                else
 
431
                  im1 = il1
 
432
                endif
361
433
                status = bas_cn2bfr(basis,ish2,if2,il2)
362
434
                if2 = if2 - iab2f + 1
363
435
                il2 = il2 - iab2f + 1
364
 
 
365
 
                nint = ( il1 - if1 + 1 ) * ( il2 - if2 + 1 )
366
 
 
367
 
                ic=1
368
 
                do iat3 = 1, nat
369
 
                  do icart = 1, 3
370
 
                    do ip1 = if1, il1
371
 
                      do ip2 = if2, il2
372
 
                        H(ic)=0.D0
373
 
                        ic = ic + 1
374
 
                      enddo
375
 
                    enddo
376
 
                  enddo
 
436
c
 
437
c               Same Cartesian vs spherical harmonic catastrophy as
 
438
c               for ish1.
 
439
c
 
440
                status = bas_continfo(basis,ish2,ityp2,nprim,ngen,
 
441
     +                                sphcart)
 
442
                if (sphcart.eq.1.and.ityp2.ge.2) then
 
443
                  im2 = if2 + (ityp2+1)*(ityp2+2)/2 - 1
 
444
                else
 
445
                  im2 = il2
 
446
                endif
 
447
 
 
448
                nint = ( im1 - if1 + 1 ) * ( im2 - if2 + 1 )
 
449
 
 
450
                do iefc = 1, nefc
 
451
 
 
452
                  ic=1
 
453
                  do iat3 = 1, 3 ! centers A, B, and C
 
454
                    do icart = 1, 3
 
455
                      do ip1 = if1, im1
 
456
                        do ip2 = if2, im2
 
457
                          H(ic)=0.0D0
 
458
                          ic = ic + 1
 
459
                        enddo
 
460
                      enddo
 
461
                    enddo
 
462
                  enddo
 
463
 
 
464
C                 1el. -cosmo- derivatives
 
465
c                 Currently calculated on for every COSMO charge
 
466
c                 separately.
 
467
                  call intd_1epot_cosmo(basis,ish1,basis,ish2,lscr,scr,
 
468
     &                 lbuf,H,dbl_mb(iefc_c+3*(iefc-1)),
 
469
     &                 dbl_mb(iefc_q+iefc-1),1)
 
470
 
 
471
C     D x H
 
472
c
 
473
c                 Do center A (associated with ish1)
 
474
c
 
475
                  ic=1
 
476
                  do icart = 1, 3
 
477
                    dE = 0.D0
 
478
                    do ip1 = if1, il1
 
479
                      do ip2 = if2, il2
 
480
                        dE = dE + dens(ip1,ip2) * H(ic)
 
481
                        ic = ic + 1
 
482
                      enddo
 
483
                    enddo
 
484
                    if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) dE = dE + dE
 
485
                    dE = dE * qfac
 
486
                    frc_cos_elq(icart,iat1)
 
487
     &              = frc_cos_elq(icart,iat1) - dE
 
488
                  enddo
 
489
c
 
490
c                 Do center B (associated with ish2)
 
491
c
 
492
                  do icart = 1, 3
 
493
                    dE = 0.D0
 
494
                    do ip1 = if1, il1
 
495
                      do ip2 = if2, il2
 
496
                        dE = dE + dens(ip1,ip2) * H(ic)
 
497
                        ic = ic + 1
 
498
                      enddo
 
499
                    enddo
 
500
                    if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) then
 
501
                      dE = dE + dE
 
502
                    else
 
503
                      dE = 0.0d0
 
504
                    endif
 
505
                    dE = dE * qfac
 
506
                    frc_cos_elq(icart,iat2)
 
507
     &              = frc_cos_elq(icart,iat2) - dE
 
508
                  enddo
 
509
c
 
510
c                 Do center C, i.e. the Cosmo charge (associated with
 
511
c                 the atom stored in int_mb(k_efciat))
 
512
c
 
513
                  do icart = 1, 3
 
514
                    dE = 0.D0
 
515
                    do ip1 = if1, il1
 
516
                      do ip2 = if2, il2
 
517
                        dE = dE + dens(ip1,ip2) * H(ic)
 
518
                        ic = ic + 1
 
519
                      enddo
 
520
                    enddo
 
521
                    if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) dE = dE + dE
 
522
                    dE = dE * qfac
 
523
                    frc_cos_elq(icart,int_mb(k_efciat+iefc-1))
 
524
     &              = frc_cos_elq(icart,int_mb(k_efciat+iefc-1)) - dE
 
525
                  enddo
 
526
 
377
527
                enddo
378
528
 
379
 
C               1el. -cosmo- derivatives
380
 
c-              call intd_1eh1(basis,ish1,basis,ish2,lscr,scr,
381
 
c-   &               lbuf,H)
382
 
 
383
 
C     D x H
384
 
 
385
 
                ic=1
386
 
                do 50, iat3 = 1, nat
387
 
                  do 40, icart = 1, 3
388
 
                    dE = 0.D0
389
 
                    do 31, ip1 = if1, il1
390
 
                      do 30, ip2 = if2, il2
391
 
                        dE = dE + dens(ip1,ip2) * H(ic)
392
 
                        ic = ic + 1
393
 
 30                   continue
394
 
 31                 continue
395
 
                    if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) dE = dE + dE
396
 
                    dE = dE * qfac
397
 
                    frc_kin(icart,iat3) = frc_kin(icart,iat3) + dE
398
 
 40               continue
399
 
 50             continue
400
 
 
401
529
 60           continue
402
530
 70         continue
403
531
 
409
537
 80     continue
410
538
 90   continue
411
539
      next = nxtask(-nproc,task_size)
 
540
c
 
541
      if (ga_nodeid().eq.0) then
 
542
c
 
543
c       Do the Nuclear - Cosmo charge part
 
544
c
 
545
        do iat1 = 1, nat
 
546
          do ich1 = 1, nefc
 
547
            dx = coords(1,iat1,geom) - dbl_mb(0+3*(ich1-1)+iefc_c)
 
548
            dy = coords(2,iat1,geom) - dbl_mb(1+3*(ich1-1)+iefc_c)
 
549
            dz = coords(3,iat1,geom) - dbl_mb(2+3*(ich1-1)+iefc_c)
 
550
            fact = -charge(iat1,geom)*dbl_mb(iefc_q+ich1-1) / 
 
551
     $             sqrt(dx*dx+dy*dy+dz*dz)**3
 
552
            dE = dx * fact
 
553
            frc_cos_nucq(1,iat1) = frc_cos_nucq(1,iat1) + dE
 
554
            frc_cos_nucq(1,int_mb(k_efciat+ich1-1)) 
 
555
     $      = frc_cos_nucq(1,int_mb(k_efciat+ich1-1)) - dE
 
556
            dE = dy * fact
 
557
            frc_cos_nucq(2,iat1) = frc_cos_nucq(2,iat1) + dE
 
558
            frc_cos_nucq(2,int_mb(k_efciat+ich1-1)) 
 
559
     $      = frc_cos_nucq(2,int_mb(k_efciat+ich1-1)) - dE
 
560
            dE = dz * fact
 
561
            frc_cos_nucq(3,iat1) = frc_cos_nucq(3,iat1) + dE
 
562
            frc_cos_nucq(3,int_mb(k_efciat+ich1-1)) 
 
563
     $      = frc_cos_nucq(3,int_mb(k_efciat+ich1-1)) - dE
 
564
          enddo
 
565
        enddo
 
566
c
 
567
c       Do cosmo charge - cosmo charge interaction
 
568
c
 
569
        invscreen = 1.0d0/(1.0d0*screen)
 
570
        do ich1 = 1, nefc
 
571
          do ich2 = 1, ich1-1
 
572
            if (ich1.ne.ich2) then
 
573
              dx = dbl_mb(0+3*(ich1-1)+iefc_c)
 
574
     $           - dbl_mb(0+3*(ich2-1)+iefc_c)
 
575
              dy = dbl_mb(1+3*(ich1-1)+iefc_c)
 
576
     $           - dbl_mb(1+3*(ich2-1)+iefc_c)
 
577
              dz = dbl_mb(2+3*(ich1-1)+iefc_c)
 
578
     $           - dbl_mb(2+3*(ich2-1)+iefc_c)
 
579
              rr = sqrt(dx*dx+dy*dy+dz*dz)
 
580
              if (rr.lt.1.0d-10) then
 
581
                fact = 0.0d0
 
582
              else
 
583
                fact = -invscreen*dbl_mb(iefc_q+ich1-1)
 
584
     $               * dbl_mb(iefc_q+ich2-1) / 
 
585
     $                 rr**3
 
586
              endif
 
587
              dE = dx * fact
 
588
              frc_cos_qq(1,int_mb(k_efciat+ich1-1)) 
 
589
     $        = frc_cos_qq(1,int_mb(k_efciat+ich1-1)) + dE
 
590
              frc_cos_qq(1,int_mb(k_efciat+ich2-1)) 
 
591
     $        = frc_cos_qq(1,int_mb(k_efciat+ich2-1)) - dE
 
592
              dE = dy * fact
 
593
              frc_cos_qq(2,int_mb(k_efciat+ich1-1)) 
 
594
     $        = frc_cos_qq(2,int_mb(k_efciat+ich1-1)) + dE
 
595
              frc_cos_qq(2,int_mb(k_efciat+ich2-1)) 
 
596
     $        = frc_cos_qq(2,int_mb(k_efciat+ich2-1)) - dE
 
597
              dE = dz * fact
 
598
              frc_cos_qq(3,int_mb(k_efciat+ich1-1)) 
 
599
     $        = frc_cos_qq(3,int_mb(k_efciat+ich1-1)) + dE
 
600
              frc_cos_qq(3,int_mb(k_efciat+ich2-1)) 
 
601
     $        = frc_cos_qq(3,int_mb(k_efciat+ich2-1)) - dE
 
602
            endif
 
603
          enddo
 
604
        enddo
 
605
      endif
 
606
c
 
607
      if (.not.ma_pop_stack(l_efciat))
 
608
     $   call errquit("grad_hnd_cos: could not deallocate l_efciat",
 
609
     $                0,MA_ERR)
 
610
      if (.not.bq_destroy(cosmo_bq_efc))
 
611
     $   call errquit("grad_hnd_cos: bq_destroy on cosmo failed",
 
612
     $                0,GEOM_ERR)
412
613
 
413
614
      return
414
615
      end