← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/nwdft/rt_tddft/propagators/exp_pseries.F

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

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

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

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/nwdft/rt_tddft/propagators/exp_pseries.F
 
1
C
 
2
C     exp_pseries.F
 
3
C
 
4
C     Brute force exponentiation of a complex operator A via a power
 
5
C     series.  Works for any square complex matrix.
 
6
C
 
7
C     e^A = \sum_{k=0}^{\infty} 1/k! A^k
 
8
C
 
9
C     We first scale A by 2^m (such that ||A|| <= 1), then square the
 
10
C     result m times to recover the answer:
 
11
C
 
12
C     [ e^{A/2^m} ]^{2^m} = e^A
 
13
C     
 
14
C     We compute the necessary m each time via the infinity norm of ||A||.
 
15
C     
 
16
      subroutine exp_pseries (params, g_za, g_zexpa)
 
17
      implicit none
 
18
 
 
19
 
 
20
#include "global.fh"
 
21
#include "errquit.fh"
 
22
#include "mafdecls.fh"
 
23
#include "stdio.fh"
 
24
#include "matutils.fh"
 
25
#include "rt_tddft.fh"
 
26
      
 
27
 
 
28
C     == Inputs ==
 
29
      type(rt_params_t), intent(in) :: params
 
30
      integer, intent(in)           :: g_za 
 
31
 
 
32
      
 
33
C     == Outputs ==
 
34
      integer, intent(in) :: g_zexpa
 
35
 
 
36
 
 
37
C     == Parameters ==
 
38
      character(*), parameter   :: pname = "exp_pseries: "
 
39
      integer, parameter        :: maxk = 1000 !max num terms in expansion
 
40
      integer, parameter        :: max_mscale = 16 !=65536, norm should never be this big
 
41
 
 
42
      
 
43
C     == Variables ==
 
44
      integer ik
 
45
      integer dtype, n1, n2
 
46
      integer g_za_scaled
 
47
      integer im
 
48
      double precision scale_val
 
49
      double complex zscale
 
50
      integer g_new_term, g_prev_term
 
51
      double complex invfac
 
52
      integer me
 
53
      logical converged
 
54
      integer num_zeroterms
 
55
      character*120 outstring
 
56
      double precision elapsed
 
57
      double precision invk
 
58
      double complex zinvk
 
59
      double precision norm
 
60
      double precision ratio
 
61
      integer mscale            !divide A by 2^m, then square result m times
 
62
 
 
63
 
 
64
      if (params%prof) call prof_start (elapsed)
 
65
 
 
66
      me = ga_nodeid ()
 
67
 
 
68
      
 
69
C
 
70
C     Check the GA.
 
71
C
 
72
      call ga_check_handle (g_za,
 
73
     $     "first argument of "//pname//" not a valid GA")
 
74
      
 
75
      call ga_inquire (g_za, dtype, n1, n2)
 
76
      
 
77
      if (dtype .ne. mt_dcpl) call errquit (
 
78
     $     pname//" only valid for complex matricies", 0, 0)
 
79
      
 
80
      if (n1 .ne. n2)
 
81
     $     call errquit (pname//"n1 must equal n2")
 
82
 
 
83
      if (.not. ga_duplicate(g_za, g_za_scaled, "g_za_scaled"))
 
84
     $     call errquit (pname//"alloc failed", 0, GA_ERR)
 
85
      if (.not. ga_duplicate(g_za, g_new_term, "g_new_term"))
 
86
     $     call errquit (pname//"alloc failed", 0, GA_ERR)
 
87
      if (.not. ga_duplicate(g_za, g_prev_term, "g_prev_term"))
 
88
     $     call errquit (pname//"alloc failed", 0, GA_ERR)
 
89
 
 
90
 
 
91
 
 
92
C
 
93
C     Determine 2^m factor and scale input matrix.
 
94
C
 
95
C     
 
96
C     We technically need to compute the scaling factor from the
 
97
C     spectral range of the W operator, but a simple norm will suffice.
 
98
C
 
99
C     If L is the largest abs eigenvalue (or the norm as we do it
 
100
C     below), then we pick scaling factor "m" such that:
 
101
C
 
102
C            2^m = 2.5*L
 
103
C     => m log 2 = log (2.5*L)
 
104
C              m = log (2.5*L) / log 2
 
105
C
 
106
C     (we multiply by 2.5 to ensure ||W|| < 1)
 
107
C
 
108
C
 
109
      call ga_norm_infinity (g_za, norm)
 
110
 
 
111
      ratio = dlog (2.5d0*norm) / dlog (2d0)
 
112
      mscale = max (int (ratio), 0)
 
113
 
 
114
 
 
115
      if (mscale .lt. 0)
 
116
     $     call errquit (pname//"negative mscale", 0, 0)
 
117
 
 
118
      if (mscale .gt. max_mscale)
 
119
     $     call errquit (pname//"mscale exceeds 16 (2^16 = 65536)",
 
120
     $     0, 0)
 
121
 
 
122
 
 
123
      scale_val = 1d0/(2d0**mscale)
 
124
      zscale = dcmplx (scale_val, 0d0)
 
125
 
 
126
      call ga_zero (g_za_scaled)
 
127
      call ga_copy (g_za, g_za_scaled)
 
128
      call ga_scale (g_za_scaled, zscale)
 
129
 
 
130
 
 
131
C
 
132
C     k = 0 term is just I
 
133
C
 
134
      call ga_zero (g_prev_term)
 
135
      call mat_set_ident (g_prev_term)
 
136
 
 
137
      call ga_zero (g_zexpa)
 
138
      call ga_copy (g_prev_term, g_zexpa)
 
139
 
 
140
      ik = 0
 
141
      converged = .false.
 
142
      num_zeroterms = 0
 
143
      
 
144
      do while (.not. converged)
 
145
         call ga_sync ()
 
146
 
 
147
         ik = ik + 1
 
148
         zinvk = z1 / dcmplx (ik)
 
149
         
 
150
         
 
151
C
 
152
C     New term is 1/k A times previous term
 
153
C
 
154
         call ga_zero (g_new_term)
 
155
         call ga_zgemm ("N", "N", n1, n1, n1,
 
156
     $        zinvk, g_prev_term, g_za_scaled, z0, g_new_term)
 
157
 
 
158
         
 
159
 
 
160
C
 
161
C     Add new term to sum and store new term as previous.
 
162
C
 
163
         call ga_add (z1, g_zexpa, z1, g_new_term, g_zexpa)
 
164
         call ga_zero (g_prev_term)
 
165
         call ga_copy (g_new_term, g_prev_term)
 
166
         
 
167
         
 
168
C
 
169
C     Check if we have converged
 
170
C
 
171
         norm = mat_norm (g_new_term)
 
172
 
 
173
         if (norm < params%tol_series) num_zeroterms = num_zeroterms + 1
 
174
 
 
175
         if (num_zeroterms .ge. params%terms_series)
 
176
     $        converged = .true.
 
177
 
 
178
 
 
179
C
 
180
C     Stop if we have failed to converge after max num terms.
 
181
C
 
182
         if (ik.gt.maxk)
 
183
     $        call errquit (pname//"failed to converge", 0, 0)
 
184
 
 
185
      enddo
 
186
 
 
187
 
 
188
C
 
189
C     Square output matrix m times to recover unscaled result.  Use
 
190
C     g_prev_term as scratch space.
 
191
C
 
192
      do im = 1, mscale
 
193
         call ga_sync ()
 
194
         call ga_zero (g_prev_term)
 
195
         call ga_zgemm ("N", "N", n1, n1, n1,
 
196
     $        z1, g_zexpa, g_zexpa, z0, g_prev_term)
 
197
         call ga_zero (g_zexpa)
 
198
         call ga_copy (g_prev_term, g_zexpa)
 
199
      enddo
 
200
 
 
201
      
 
202
C
 
203
C     Clean up
 
204
C
 
205
      if (.not. ga_destroy(g_new_term))
 
206
     $     call errquit (pname//"destroy failed", 0, GA_ERR)
 
207
      if (.not. ga_destroy(g_prev_term))
 
208
     $     call errquit (pname//"destroy failed", 0, GA_ERR)
 
209
      if (.not. ga_destroy(g_za_scaled))
 
210
     $     call errquit (pname//"destroy failed", 0, GA_ERR)
 
211
 
 
212
      
 
213
C
 
214
C     If profiling is enabled print convergence data.
 
215
C      
 
216
      if (params%prof) then 
 
217
         write(outstring,"(a,i0,a,i0,a)")
 
218
     $        "Power series with scaling 2^",
 
219
     $        mscale, " converged after ", ik," terms;"
 
220
         call prof_end (elapsed, trim(outstring))
 
221
      endif
 
222
 
 
223
      end subroutine
 
224