← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/nwdft/rt_tddft/init/rt_tddft_init_geoms.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/init/rt_tddft_init_geoms.F
 
1
C
 
2
C     Get geometry info from rtdb and load into struct.
 
3
C
 
4
C     Each geometry (fragment) has a "map" which is nbf_ao x nbf_ao, where nbf_ao
 
5
C     is the number of basis functions in the "active" geometry (ie set
 
6
C     geometry in input deck).  You can then project the full nbf_ao x nbf_ao
 
7
C     dipole integrals, etc onto each individual fragment via
 
8
C     ga_elem_multiply ().
 
9
C
 
10
C     NOTE: params%geom_active_handle refers to the nwchem geom
 
11
C     subsystem handle for "geometry", while params%active_geom refers
 
12
C     to the params%geom(:) array of geom structs.
 
13
C
 
14
      subroutine rt_tddft_init_geoms (params)
 
15
      implicit none
 
16
 
 
17
#include "errquit.fh"
 
18
#include "mafdecls.fh"
 
19
#include "stdio.fh"
 
20
#include "global.fh"
 
21
#include "msgids.fh"
 
22
#include "util.fh"
 
23
#include "rtdb.fh"
 
24
#include "cdft.fh"
 
25
#include "inp.fh"
 
26
#include "geomP.fh"
 
27
#include "geom.fh"
 
28
#include "bas.fh"
 
29
#include "sym.fh"
 
30
#include "rt_tddft.fh"
 
31
 
 
32
C
 
33
C     Note: geomP.fh supplies ngeom_rtdb, names_rtdb, max_geom_rtdb, ndipole
 
34
C
 
35
      
 
36
 
 
37
C     == In/out ==
 
38
      type(rt_params_t) params  !geom params stored in here
 
39
 
 
40
 
 
41
C     == Parameters ==
 
42
      character(*), parameter :: pname = "rt_tddft_init_geoms: "
 
43
      character(*), parameter :: geom_tmp_name = "rt_tmp_geom"
 
44
      double precision, parameter :: dist_thresh = 1d-8
 
45
 
 
46
      
 
47
C     == Variables ==
 
48
      integer me
 
49
      integer ig
 
50
      
 
51
      integer geom_tmp
 
52
      character*16 geom_rtdb_name, geom_excite_name
 
53
      character*16 active_geom_name
 
54
      integer geom_tmp_ncent
 
55
 
 
56
      integer ibf, icen, jcen
 
57
      character*16 icen_tag, jcen_tag
 
58
      double precision icen_loc(3), jcen_loc(3)
 
59
      double precision icen_charge, jcen_charge
 
60
      double precision dist
 
61
      logical bf_in_geom
 
62
      integer nbf_geom
 
63
      integer per_bf, per_cen
 
64
      double precision geom_charge_nuc
 
65
      logical got_excite
 
66
      integer is
 
67
 
 
68
C      character*25 role_string
 
69
      character*2 role_string
 
70
 
 
71
      logical found_active
 
72
 
 
73
      double precision ndip(3)
 
74
 
 
75
 
 
76
      me = ga_nodeid ()
 
77
 
 
78
 
 
79
C
 
80
C     Check that rt hardcoded max num geoms is less than = nwchem hardlimit
 
81
C
 
82
      if (rt_max_geoms .lt. nw_max_geoms) !hardcoded in nwc_const.fh
 
83
     $     call errquit (pname//"rt_max_geoms < nw_max_geoms",0,0)
 
84
 
 
85
 
 
86
C
 
87
C     Get geom info from rtdb
 
88
C
 
89
      if (.not.rtdb_get (params%rtdb, "geometry:ngeom",
 
90
     $     mt_int, 1, ngeom_rtdb))
 
91
     $     call errquit (pname//"failed to get ngeom from rtdb",0,0)
 
92
      
 
93
      if (ngeom_rtdb .lt. 1)
 
94
     $     call errquit (pname//"invalid ngeom_rtdb < 1",0,0)
 
95
 
 
96
c$$$      if (ngeom_rtdb .gt. nw_max_geoms)
 
97
c$$$     $     call errquit (pname//"invalid ngeom_rtdb > nw_max_geoms",0,0)
 
98
 
 
99
      if (ngeom_rtdb .gt. max_geom_rtdb)
 
100
     $     call errquit(pname//"invalid ngeom_rtdb > max_geom_rtdb",0,0)
 
101
 
 
102
      if (.not. rtdb_cget (params%rtdb,'geometry:names',
 
103
     $     ngeom_rtdb,names_rtdb))
 
104
     $     call errquit (pname//"failed to read names from rtdb",0,0)
 
105
 
 
106
 
 
107
C
 
108
C     Get internal index for active geometry from the rtdb name for the
 
109
C     geometry.
 
110
C
 
111
 
 
112
      if (.not. rtdb_cget (params%rtdb,'geometry',1,active_geom_name))
 
113
     $     call errquit (pname//"failed to read active geom from rtdb",
 
114
     $     0,0)
 
115
 
 
116
      found_active = .false.
 
117
      
 
118
      do ig = 1, ngeom_rtdb
 
119
         if ( trim(names_rtdb(ig)) .eq. trim(active_geom_name) ) then
 
120
            if (found_active) then
 
121
               call errquit (pname//"found multiple active geoms",0,0)
 
122
            else
 
123
               params%geom_active = ig
 
124
               found_active = .true.
 
125
            endif
 
126
         endif
 
127
      enddo
 
128
 
 
129
      if (.not. found_active)
 
130
     $     call errquit (pname//"failed to find active geom",0,0)
 
131
 
 
132
 
 
133
C
 
134
C     Loop over all geoms found, create a temporary geometry object for
 
135
C     each, calculate some properties, store in the rt struct, and
 
136
C     delete the geom object (leaving geoms untouched).
 
137
C
 
138
      do ig = 1, ngeom_rtdb
 
139
 
 
140
C     (get from rtdb)
 
141
         if (.not.geom_create (geom_tmp, geom_tmp_name))
 
142
     $        call errquit (pname//"geom_create failed",0,0)
 
143
         
 
144
         if (.not. geom_check_handle (geom_tmp,pname))
 
145
     $        call errquit (pname//"not a valid geom",0,0)
 
146
 
 
147
         if (.not. geom_rtdb_load (params%rtdb,
 
148
     $        geom_tmp, names_rtdb(ig)))
 
149
     $        call errquit (pname//"failed to load geom",0,0)
 
150
 
 
151
C     (extract info from geom)
 
152
         if (.not. geom_ncent (geom_tmp, geom_tmp_ncent))
 
153
     $        call errquit (pname//"geom_ncent failed",0,0)
 
154
 
 
155
         if (.not. geom_nuc_dipole (geom_tmp, ndip))
 
156
     $        call errquit (pname//"geom_nucl_dipole failed",0,0)
 
157
 
 
158
         geom_rtdb_name = names_rtdb(ig)
 
159
 
 
160
C     
 
161
C     Check that geometry fragment is not larger than active geometry.
 
162
C
 
163
         if (geom_tmp_ncent .gt. params%natoms) call errquit (pname//
 
164
     $        'geom "'//trim(geom_rtdb_name)//
 
165
 
 
166
     $        '" has more atoms than active geometry',
 
167
     $        geom_tmp_ncent, GEOM_ERR)
 
168
 
 
169
         
 
170
C     (alloc and compute mask GA)
 
171
         if (.not.ga_create (mt_dbl, params%nbf_ao, params%nbf_ao,
 
172
     $        trim(geom_rtdb_name)//"_mask", 0, 0,
 
173
     $        params%geom(ig)%g_mask))
 
174
     $        call errquit (pname//"alloc mask failed", 0, GA_ERR)
 
175
 
 
176
 
 
177
C
 
178
C     Mask starts as all 1's, then we 0 any row/column which refers to a
 
179
C     basis function not in this geom.
 
180
C
 
181
         call ga_fill (params%geom(ig)%g_mask, 1d0)
 
182
         
 
183
 
 
184
C
 
185
C     Loop over basis functions in overall/active basis set.  For each
 
186
C     basis function, find the "owning" atom, and then check if that
 
187
C     atom is part of the "ig" geometry.
 
188
C
 
189
         nbf_geom = 0
 
190
         do ibf = 1, params%nbf_ao
 
191
            if (.not. bas_bf2ce (params%ao_bas_han, ibf, icen))
 
192
     $           call errquit (pname//"bas_bf2ce failed", 0, 0)
 
193
            
 
194
C     (note this acts on full active geom, specified by the handle
 
195
C     stored in params)
 
196
            if (.not. geom_cent_get (params%geom_active_handle, icen,
 
197
     $           icen_tag, icen_loc, icen_charge))
 
198
     $           call errquit (pname//"geom_cent_get active failed",0,0)
 
199
 
 
200
            bf_in_geom = .false.
 
201
            do jcen = 1, geom_tmp_ncent  !increment over atoms in "ig" geom
 
202
 
 
203
               if (.not. geom_cent_get (geom_tmp, jcen,
 
204
     $              jcen_tag, jcen_loc, jcen_charge))
 
205
     $              call errquit(pname//"geom_cent_get jcen failed",0,0)
 
206
               
 
207
C     (icen_loc = location of ibf; jcen_loc = location of jcen'th atom/center in geom)
 
208
               dist = dsqrt (
 
209
     $              ( icen_loc(1) - jcen_loc(1) )**2 +
 
210
     $              ( icen_loc(2) - jcen_loc(2) )**2 +
 
211
     $              ( icen_loc(3) - jcen_loc(3) )**2 )
 
212
 
 
213
               if (dist < dist_thresh) bf_in_geom = .true.
 
214
               
 
215
            enddo ! jcen loop
 
216
 
 
217
 
 
218
C     (this basis func not in this geom, zero corresponding rows+columns)
 
219
            if (.not. bf_in_geom) then
 
220
               call ga_fill_patch (params%geom(ig)%g_mask,
 
221
     $              1, params%nbf_ao, ibf, ibf, 0d0)
 
222
               call ga_fill_patch (params%geom(ig)%g_mask,
 
223
     $              ibf, ibf, 1, params%nbf_ao, 0d0)
 
224
            else ! this bf in this geom
 
225
               nbf_geom = nbf_geom + 1
 
226
            endif
 
227
         enddo ! ibf loop
 
228
 
 
229
         
 
230
C     (compute total nuclear charge on this geom)
 
231
         geom_charge_nuc = 0d0
 
232
         do jcen = 1, geom_tmp_ncent
 
233
 
 
234
            if (.not. geom_cent_get (geom_tmp, jcen,
 
235
     $           jcen_tag, jcen_loc, jcen_charge))
 
236
     $           call errquit(pname//"geom_cent_get jcen failed",0,0)
 
237
            
 
238
            geom_charge_nuc =  geom_charge_nuc + jcen_charge
 
239
         enddo
 
240
 
 
241
 
 
242
C     (store misc params)
 
243
         params%ngeoms = ngeom_rtdb
 
244
         params%geom(ig)%name  = trim (geom_rtdb_name)
 
245
         params%geom(ig)%ncent = geom_tmp_ncent
 
246
         params%geom(ig)%nbf = nbf_geom
 
247
         params%geom(ig)%charge_nuc = geom_charge_nuc
 
248
         params%geom(ig)%ndip%x = ndip(1)
 
249
         params%geom(ig)%ndip%y = ndip(2)
 
250
         params%geom(ig)%ndip%z = ndip(3)
 
251
         params%geom(ig)%measure = .true. !XXX HARDCODED MEASURE ALL
 
252
 
 
253
         
 
254
C     (remove tmp geom object, rtdb untouched)
 
255
         if (.not. geom_destroy (geom_tmp))
 
256
     $        call errquit (pname//"geom_destroy failed", 0, 0)
 
257
 
 
258
      enddo ! ig loop
 
259
 
 
260
      
 
261
C
 
262
C     Print geom info to stdout
 
263
C
 
264
      if (me .eq.0) then
 
265
         write (luout, *) " "
 
266
         write (luout, *) " "
 
267
         
 
268
         write (luout, *) "Geometry               "//
 
269
     $        "Atoms      Basis func.   "//
 
270
     $        "Nuc. charge          Nuc. dip. mom."
 
271
 
 
272
         write (luout, *) "-----------------"//
 
273
     $        "------------------------------------"//
 
274
     $        "------------------------------------"
 
275
 
 
276
      endif
 
277
 
 
278
      do ig = 1, params%ngeoms
 
279
 
 
280
         per_bf = 100 * params%geom(ig)%nbf / params%nbf_ao
 
281
         per_cen = 100 * params%geom(ig)%ncent / params%natoms
 
282
 
 
283
C         role_string = " "
 
284
 
 
285
C         if ( trim(params%geom(ig)%name) .eq. trim(active_geom_name))
 
286
C     $        role_string = trim(role_string) // "<= Active "
 
287
 
 
288
         if (ig .eq. params%geom_active) then
 
289
            role_string = " *"
 
290
         else
 
291
            role_string = "  "
 
292
         endif
 
293
 
 
294
         if (me.eq.0) then
 
295
 
 
296
            write (luout,
 
297
C     $           "(2x,a16,2x,i6,a,i3,a,5x,i6,a,i3,a,4x,1f8.2,4x,a)")
 
298
     $           "(a,a16,i6,a,i3,a,i6,a,i3,a,2x,1f8.2,2x,3f10.2)")
 
299
     $           role_string,
 
300
     $           params%geom(ig)%name,
 
301
     $           params%geom(ig)%ncent, " (", per_cen , "% )", 
 
302
     $           params%geom(ig)%nbf,   " (", per_bf, "% )",
 
303
     $           params%geom(ig)%charge_nuc,
 
304
     $           params%geom(ig)%ndip%x,
 
305
     $           params%geom(ig)%ndip%y,
 
306
     $           params%geom(ig)%ndip%z
 
307
         endif
 
308
         
 
309
         call util_flush (luout)
 
310
      enddo
 
311
 
 
312
 
 
313
C
 
314
C     Now that we have determined the active geometry, update any
 
315
C     excitations with -1 geom (ie if exciting orbitals directly) to
 
316
C     excite the entire active geometry.
 
317
C
 
318
c$$$      do is = 1, params%nexcites
 
319
c$$$         if (params%excite(is)%geom_indx .lt. 1) then
 
320
c$$$            params%excite(is)%geom_indx = params%geom_active
 
321
c$$$         endif
 
322
c$$$      enddo
 
323
 
 
324
 
 
325
 
 
326
C
 
327
C     It appears that dipole integrals only work for Abelian symmetry,
 
328
C     so we enforce that for each geom.
 
329
C
 
330
C
 
331
      if (.not.sym_abelian_group(params%geom_active_handle))
 
332
     $     call errquit (pname//
 
333
     $     "active geom: non-Abelian symmetry; disable autosym",
 
334
     $     0, GEOM_ERR)
 
335
 
 
336
C     (debug)
 
337
C      call util_flush (luout)
 
338
C      if (.not. bas_print_all())
 
339
C     $     call errquit (pname//"basis print failed")
 
340
C      call util_flush (luout)
 
341
      
 
342
      end subroutine
 
343