← Back to branch summary

~ubuntu-branches/ubuntu/saucy/nwchem/saucy

« back to all changes in this revision

Viewing changes to .pc/09_backported_6.1.1_fixes.patch/src/property/raman.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2012-09-13 14:35:37 UTC
  • Revision ID: package-import@ubuntu.com-20120913143537-5gkt1l68rxrprgnl
Tags: 6.1-4
[ Michael Banck ]
* debian/patches/09_backported_6.1.1_fixes.patch: New patch, backports the
  source code changes of the nwchem-6.1.1 bugfix release:
  + PW: Fixed backspace issues on file I/O that caused I/O errors.
  + DFT: Removed dummy and bq centers from the Grimme dispersion
    corrections.
  + DFT: Fixed a race condition in the density fitting code.
  + DFT: Added a check for singularities in the HCTH functionals.
  + DFT: Fixed a problem with the DFT grids which caused strange behaviors
    if the number of cores is so large that some cores do not get any grid
    points.
  + HF&DFT: Fixed rolling back to distributed memory Fock-builder if not
    enough memory is available to use the replicated data one. Previously
    the code would crash trying to use non-existing GAs.
  + HF&DFT: Fixed clashes between MPI and GA communication when using OpenIB
    which enhances the stability.
  + MP2&DFT: On systems with limited I/O capabilities some quantities like
    2-electron integrals and DFT grids are now stored in memory rather than
    on disk.
  + CASSCF: Added ga_sync to fix race conditions that can cause the Davidson
    diagonalizer to fail.
  + CASSCF: Fixed a problem with the phase in the Lagrangian that caused
    problems with the gradient evaluation.
  + RAMAN: A number of problems with static polarizabilities were fixed.
  + Property: Fixed an issue with add_patch that caused unexpected results
    with dynamic polarizabilities.
  + DRDY: Removed system calls to copy files avoiding forking from NWChem
    processes which is relatively likely to fail due to the resources
    attached to such a process.
  + Input: Fixed some issues with GEOM LOAD that caused the selection of
    centers to fail in some cases.
  + Geometry: Dummy centers are no longer removed from a geometry so that
    constraints involving those centers remain valid.
  + Memory: All shared memory (global memory region) is now allocated at the
    start.

[ Daniel Leidert ]
* debian/control: Added X-Python-Version.
  (Build-Depends): Added python-dev. Use texlive to fix manual build.
  (Standards-Version): Bumped to recent 3.9.3.
* debian/nwchem.1: Added.
* debian/nwchem.doc-base: Ditto.
* debian/nwchem.lintian-overrides: Ditto.
* debian/nwchem.manpages: Ditto.
* debian/nwchem-data.lintian-overrides: Ditto.
* debian/rules: Added PYTHONVERSION, PYTHONHOME. Enable parallel building.
  (NWCHEM_MODULES): Build with python support (pnnl).
* debian/patches/02_makefile_flags.patch: Adjusted.
  - src/config/makefile.h: Fix linker flags building with python support.

Show diffs side-by-side

added added

removed removed

Lines of Context:
.pc/09_backported_6.1.1_fixes.patch/src/property/raman.F
 
1
c
 
2
c ---------------------------------------------------------------------
 
3
c computes the xyz coords along the imode normal mode
 
4
c ---------------------------------------------------------------------
 
5
c
 
6
      subroutine raman_modestep(rtdb,nat3,natom,geom,rmmodes,imode,
 
7
     &    iii,first,eigenvecs,eigenvals,ncoords,rminfo,step_size)
 
8
c
 
9
      implicit none
 
10
c
 
11
#include "errquit.fh"
 
12
#include "rtdb.fh"
 
13
#include "mafdecls.fh"
 
14
#include "stdio.fh"
 
15
#include "geom.fh"
 
16
#include "nwc_const.fh"
 
17
#include "inp.fh"
 
18
#include "global.fh"
 
19
#include "bas.fh"
 
20
c
 
21
      integer iii ! step number which determines sign 
 
22
      integer imode ! mode cordinates are dipsplaced along
 
23
      integer iatom, ixyz, ivec ! counting index
 
24
      integer geom    ! [input] geom handle
 
25
      integer rtdb    ! [input] rtdb handle
 
26
      integer natom   ! [input] number of atoms
 
27
      integer nat3    ! [input] 3*number of atoms
 
28
      integer first   ! first mode to consider in aoresponse (default =7 ramana =1 hyperraman)
 
29
      integer tmpmode ! set to fill rminfo from 1 ( not 7 for raman calc)
 
30
      integer rmmodes ! # of raman active modes
 
31
c
 
32
      double precision rminfo(rmmodes,4) ! data for raman spec
 
33
      double precision step_size ! multiplictive factor for step along normal mode
 
34
      double precision sign,bohr2ang ! sign of the step , unit conversion
 
35
      double precision eigenvecs(nat3,nat3) ! [input](xyz&atom,mode)
 
36
      double precision eigenvals(nat3)      ! [input] (mode)
 
37
      double precision ncoords(3,natom)    ! [scratch] coords after step
 
38
      double precision steps(3,natom)     ! [scratch] step generated by vector and scaled
 
39
c
 
40
      parameter (bohr2ang=0.52917724924D+00) ! CONVERSION OF BOHR to ANGSTROMS
 
41
c -------------determine sign of the step---------------------------------
 
42
      if (iii.eq.1) then
 
43
        sign = 1.0D+00
 
44
      else
 
45
        sign = -1.0D+00
 
46
      endif
 
47
c  --- save imode values in rminfo  ---
 
48
      tmpmode=imode-first+1
 
49
      rminfo(tmpmode,1) = eigenvals(imode)
 
50
c --------------------------------------------------------------------
 
51
        ivec = 1
 
52
        do iatom = 1,natom
 
53
          do ixyz = 1,3
 
54
      steps(ixyz,iatom)=sign*step_size*eigenvecs(ivec,imode)
 
55
            ivec = ivec + 1
 
56
          enddo ! ixyz
 
57
        enddo ! iatom
 
58
c
 
59
      call daxpy(nat3,1.0d00,steps,1,ncoords,1) ! mult coords 
 
60
      if (.not. geom_cart_coords_set(geom,ncoords))
 
61
     $     call errquit('raman_modestep: bad geom',0, GEOM_ERR)
 
62
c
 
63
      return
 
64
      end
 
65
c
 
66
c --------------------------------------------------------------------
 
67
c    Finite differencing for Raman Scattering Calculations
 
68
c --------------------------------------------------------------------
 
69
c
 
70
      subroutine fd_raman(rtdb,imode,rmmodes,natom,nat3,alfarex2,
 
71
     &           alfarimx2,step_size,rminfo,eigenvecs,mass)
 
72
c
 
73
      implicit none
 
74
c
 
75
#include "errquit.fh"
 
76
#include "global.fh"
 
77
#include "mafdecls.fh"
 
78
#include "stdio.fh"
 
79
#include "rtdb.fh"
 
80
#include "util.fh"
 
81
c
 
82
      logical debug
 
83
c
 
84
      integer rmmodes !  # modes used in raman calculation (3N-6)
 
85
      integer  i,j,m,ii,jj,mm,ivec  ! counting indexes
 
86
      integer rtdb    ! [input] rtdb handle
 
87
      integer imode   ! mode #
 
88
      integer natom   ! [input] number of atoms
 
89
      integer nat3    ! [input] 3*number of atoms
 
90
c
 
91
      double precision rminfo(rmmodes,4) ! raman data
 
92
      double precision step_size,stepsize ! [input] step of finite differencing
 
93
c      double precision laser(nfreq) ! [input] frequency used in AORESPONSE 
 
94
      double precision eigenvecs(nat3,nat3) ! [input](xyz&atom,mode)
 
95
      double precision tmode(3,natom) ! [input](atom#,xyz)
 
96
      double precision mass(natom) ! [input](atom#)
 
97
      double precision norm ! constant
 
98
      double precision fdipol(3,3), fdrpol(3,3) ! finite difference 
 
99
      double precision alfarex2(6,3), alfarimx2(6,3) ! AORESPONSE date for plus and minus step
 
100
      double precision ar2, ai2, gr2, gi2 ! alpha and g for imode
 
101
      double precision zero, one, two, three, six, seven,
 
102
     &                  fourty_five, bohr2ang
 
103
      parameter (zero=0.0D+00, one=1.0D+00, two=2.0D+00, three=3.0D+00,
 
104
     &            six=6.0D+00, seven=7.0D+00, fourty_five=45.0D+00)
 
105
      parameter (bohr2ang=0.52917724924D+00) ! CONVERSION OF BOHR to ANGSTROMS
 
106
c
 
107
      debug =  ( .false. .and. ga_nodeid().eq.0 )
 
108
c --------------------------------------------------------------------  
 
109
      call dfill(9,0.0D+00,fdrpol,1)  ! real fdipole polariability for 
 
110
      call dfill(9,0.0D+00,fdipol,1) ! imaginary fdipole polariability for
 
111
      call dfill(3*natom,0.0D+00,tmode,1) ! 
 
112
c zero
 
113
      stepsize = zero
 
114
      m = imode - 6 
 
115
      j=1
 
116
      i=1
 
117
      ar2   = zero ! alpha real
 
118
      gr2   = zero ! gradient real
 
119
      ai2   = zero ! alpha imaginaruy
 
120
      gi2   = zero ! gradient imaginary
 
121
      MM=3
 
122
c -----------mass weight the normal coordinates--------------------------------
 
123
       norm= 0.0D+00
 
124
       ivec = 1
 
125
      do i=1,natom
 
126
       do j=1,3
 
127
         tmode(j,i) = eigenvecs(ivec,imode)*sqrt(mass(i))
 
128
          ivec = ivec + 1
 
129
       end do
 
130
      end do
 
131
c ---calculate the norm of the qmass---
 
132
      do i=1,natom
 
133
         norm = norm + tmode(1,i)*tmode(1,i) + 
 
134
     *    tmode(2,i)*tmode(2,i) + tmode(3,i)*tmode(3,i)
 
135
      end do
 
136
      norm = sqrt(norm)
 
137
c --- normalize ---
 
138
      do i=1,natom
 
139
       do j=1,3
 
140
         tmode(j,i)= tmode(j,i)/norm
 
141
       end do
 
142
      end do
 
143
c
 
144
      do i=1,natom
 
145
       do j=1,3
 
146
         tmode(j,i) = tmode(j,i)/sqrt(mass(i))
 
147
       end do
 
148
      end do
 
149
c
 
150
      norm = zero
 
151
      do i=1,natom
 
152
         norm = norm + tmode(1,i)*tmode(1,i) + tmode(2,i)*tmode(2,i)
 
153
     *          + tmode(3,i)*tmode(3,i)
 
154
      end do
 
155
c --- calculate stepsize --- 
 
156
      stepsize = step_size/sqrt(norm)
 
157
c
 
158
      if (ga_nodeid().eq.0)  write(luout,*) 'Stepsize :',stepsize
 
159
c --------------------------------------------------------------------
 
160
      j=1
 
161
      i=1
 
162
      DO ii=1,3
 
163
       DO jj=1,3
 
164
c       difference polaraizability (real)
 
165
      fdrpol(ii,jj)     = alfarex2(ii,jj) - alfarex2(ii+3,jj)
 
166
c        difference polaraizability (imaginary)
 
167
      fdipol(ii,jj)     = alfarimx2(ii,jj) - alfarimx2(ii+3,jj)
 
168
c        convert units
 
169
       fdrpol(ii,jj) = ( fdrpol(ii,jj)*(bohr2ang**2) )/(two*stepsize)
 
170
       fdipol(ii,jj) = ( fdipol(ii,jj)*(bohr2ang**2) )/(two*stepsize)
 
171
       enddo
 
172
      enddo
 
173
c ----------calculate a and g for real and impaginary--------
 
174
      ar2 = ( ( (fdrpol(i,j)+fdrpol(i+1,j+1)+fdrpol(i+2,j+2))/three) 
 
175
     *       *( (fdrpol(i,j)+fdrpol(i+1,j+1)+fdrpol(i+2,j+2))/three) ) 
 
176
c
 
177
      gr2 = ( ( ((fdrpol(i,j)    - fdrpol(i+1,j+1))**2)
 
178
     *      +  ((fdrpol(i+1,j+1) - fdrpol(i+2,j+2))**2) 
 
179
     *      +  ((fdrpol(i+2,j+2) - fdrpol(i,j    ))**2) )
 
180
     *      + ( six*((fdrpol(i,j+1)   * fdrpol(i,j+1))
 
181
     *      +        (fdrpol(i+2,j)   * fdrpol(i+2,j))
 
182
     *      +        (fdrpol(i+1,j+2) * fdrpol(i+1,j+2))) )  )/two
 
183
c imaginary
 
184
      ai2 = ( ( (fdipol(i,j)+fdipol(i+1,j+1)+fdipol(i+2,j+2))/three)   
 
185
     *       *( (fdipol(i,j)+fdipol(i+1,j+1)+fdipol(i+2,j+2))/three) ) 
 
186
c
 
187
      gi2 = (  ( ((fdipol(i,j)     - fdipol(i+1,j+1))**2)
 
188
     *      +    ((fdipol(i+1,j+1) - fdipol(i+2,j+2))**2)
 
189
     *      +    ((fdipol(i+2,j+2) - fdipol(i,j    ))**2) )
 
190
     *      + ( six*((fdipol(i,j+1)   * fdipol(i,j+1))
 
191
     *      +        (fdipol(i+2,j)   * fdipol(i+2,j))
 
192
     *      +        (fdipol(i+1,j+2) * fdipol(i+1,j+2))) )  )/two
 
193
c --------------------------------------------------------------------
 
194
      if (debug) then
 
195
          write(luout,*) 'alfarimx2'
 
196
          call output(alfarimx2,1,6,1,3,6,3,1)
 
197
          write(luout,*) 'alfarex2'
 
198
          call output(alfarex2,1,6,1,3,6,3,1)
 
199
          write(luout,*) 'fdrpol'
 
200
          call output(fdrpol,1,3,1,3,3,3,1)
 
201
          write(luout,*) 'fdipol'
 
202
          call output(fdipol,1,3,1,3,3,3,1)
 
203
          write(luout,*) 'ar2=',ar2
 
204
          write(luout,*) 'gr2=',gr2
 
205
          write(luout,*) 'ai2=',ai2
 
206
          write(luout,*) 'gi2=',gi2
 
207
          write(luout,*) 'stepsize=',stepsize
 
208
      endif !  for DEBUG
 
209
      rminfo(m,2) = fourty_five*ar2 + seven*gr2 ! real component  
 
210
      rminfo(m,3) = fourty_five*ai2 + seven*gi2 ! imaginary component
 
211
      rminfo(m,4) = fourty_five*ar2 + seven*gr2 + fourty_five*ai2 +
 
212
     &               seven*gi2 ! total fdipole_polarizability_derivative
 
213
c------------ print result from finite difference ---------------------
 
214
      if (ga_nodeid().eq.0) write(luout,9001) rminfo(m,1),rminfo(m,2)
 
215
      if (ga_nodeid().eq.0) write(luout,9002) rminfo(m,1),rminfo(m,3)
 
216
      if (ga_nodeid().eq.0) write(luout,9003) rminfo(m,1),rminfo(m,4)
 
217
c --------------- save the rminfo for restart -------------------------
 
218
      if (.not. rtdb_put(rtdb, 'raman:rminfo ', mt_dbl,
 
219
     &    rmmodes*4, rminfo))
 
220
     &  call errquit('FD_raman:failed to write rminfo', 0, RTDB_ERR)
 
221
c --------------------------------------------------------------------
 
222
 9001 FORMAT(/1X,F7.2,2X,'REAL     ',2X,F24.6/)
 
223
 9002 FORMAT(/1X,F7.2,2X,'IMAGINARY',2X,F24.6/)
 
224
 9003 FORMAT(/1X,F7.2,2X,'TOTAL    ',2X,F24.6/)
 
225
      return
 
226
      end 
 
227
c
 
228
c --------------------------------------------------------------------
 
229
c   calculate absolute differential raman scattering cross section.
 
230
c   Producing a raman scattering output file 
 
231
c --------------------------------------------------------------------
 
232
c
 
233
      subroutine raman_scattering(rtdb,start,last,rmmodes,nfreq,plot,
 
234
     &                            line,width,laser,rminfo)
 
235
c
 
236
      implicit none
 
237
c
 
238
#include "global.fh"
 
239
#include "stdio.fh"
 
240
#include "rtdb.fh"
 
241
#include "inp.fh"
 
242
c
 
243
      character*16 plot, line
 
244
c
 
245
      integer rmmodes ! number or raman active modes
 
246
      integer rtdb    ! [input] rtdb handle
 
247
      integer nfreq   ! [input] number of excitation freqencies
 
248
      integer i,j   ! counting indexes 
 
249
      integer numpts ! number of points in plot
 
250
      integer last  ! [input] last mode to calculate for aoresponse
 
251
      integer begin ! first=begin if not restarting, else modified
 
252
      integer start ! [input] first mode to use in calculation of plot
 
253
c
 
254
      integer ifreq
 
255
      Double precision laser(nfreq) ! [input] laser used in AORESPONSE calculation
 
256
      Double precision width ! [input] line width
 
257
      Double precision rminfo(rmmodes,4) ! raman data
 
258
c      Double precision freq(rmmodes) ! [input] number of excitation frequencies used in AORESPONSE
 
259
      Double precision tempfc, frq4th 
 
260
      Double precision crs_real(rmmodes) ! real differential cross section
 
261
      Double precision  crs_imag(rmmodes) ! imaginary differential cross section
 
262
      Double precision crs_tot(rmmodes) ! total differential cross section
 
263
      Double precision frq, frqmin, frqmax, dfrq ! current, minimum, maximum and step of frequency
 
264
      Double precision afac, bfac, factor ! lineshape factors.
 
265
      Double precision conv_tot, conv_real, conv_imag !
 
266
      Double precision lambda ! exciation inverse wavenumber (1/cm)
 
267
      Double precision widthfac, pi, planck, boltz, speed, epsilon0 
 
268
      Double precision avogadro, nm2icm, temp, scale, amu,au2nm
 
269
      Double precision exparg, conver, zero, two, one, forty_five
 
270
c
 
271
      character*255 filename
 
272
      integer unitno
 
273
      parameter (unitno = 77)
 
274
c
 
275
      parameter (widthfac = 2.35482D+00)
 
276
      parameter (zero = 0.0D+00, one = 1.0D+00, 
 
277
     *           two = 2.0D+00, forty_five = 45.0D+00)
 
278
      parameter (pi = 3.14159265358979323846D+00)
 
279
      parameter (planck = 6.6260755D-34) ! plank's constant (J*sec)
 
280
      parameter (boltz = 1.3806580D-23) ! bolztmann's constant  (J/K)
 
281
      parameter (speed = 2.99792458D+08) ! speed of light (m/sec)
 
282
      parameter (epsilon0 = 8.8541878D-12) ! free permitivity (J^-1C^2m^-1)
 
283
      parameter (avogadro = 6.02214199D+23) ! avagadro's number (unitless)
 
284
      parameter (amu = (one/avogadro)*1.0D-03 ) ! atomic mass unit 
 
285
      parameter (au2nm =45.563353D+00)  ! converion factor from au to nm
 
286
      parameter (nm2icm = 1.0D+07) !converion factor from nm to 1/cm
 
287
      parameter (temp = 3.0D+02) ! tempature (k)
 
288
      parameter (scale = 1.0D+34)
 
289
      parameter (exparg  = (planck*speed*100.0D+00)/boltz )
 
290
      parameter (conver  = ( (two*planck*pi*pi)/speed*1.0D-40)/amu ) ! conver = 2*planck_h*pi^2/c * 10^40/atomic_mass_unit               
 
291
      parameter (numpts  = 1000)
 
292
c
 
293
      logical debug
 
294
      debug      = .false. ! produce debug output
 
295
c
 
296
c ----------outline and units of the calculation:---------------------
 
297
c d(sigma)/d(omega) = h/(8*(epsilon0)^2 *c) *(lambda-lambda_z)^4/lambda_z
 
298
c                     /(1-exp(-h*c*lambda_z/(kT))*(scattering factor)/45                    
 
299
c differential cross sections:
 
300
c - constant pre-factor for differential scattering cross sections:
 
301
c        h/(8*(epsilon0)^2*c) = 3.524100053d-21 J^3*s^2*m/C^4
 
302
c - lambda_t is given in [1/m^3]
 
303
c - squared polarizability derivatives have to be given in
 
304
c   [C^2*m^2/V^2/kg], but they are given as squared derivatives of
 
305
c   polarizability volumes in [Angstrom^4/amu].
 
306
c   conversion factor: (4*pi*epsilon0)^2*10^-40{Angstrom^4/m^4}
 
307
c                      /atomic_mass_unit{amu/kg}
 
308
c --------------------------------------------------------------------
 
309
c
 
310
      if (debug) then
 
311
          write(6,*) "rmmodes",rmmodes
 
312
          do ifreq = 1,nfreq
 
313
            write(6,*) "laser",laser(ifreq)
 
314
          end do
 
315
          write(6,*) 'rminfo'
 
316
          call output(rminfo,1,rmmodes,1,4,rmmodes,4,1)
 
317
      endif
 
318
c
 
319
c     --- write raman scattering data to file ---
 
320
c
 
321
      call util_file_name(plot,.false.,.false.,filename)
 
322
      open(unitno, status='unknown', form='formatted',file=filename)
 
323
 
 
324
      write(unitno,9002) plot
 
325
      do ifreq = 1,nfreq
 
326
         write(unitno,9003) au2nm/laser(ifreq)
 
327
      end do
 
328
      write(unitno,9004) line 
 
329
      write(unitno,9005) width
 
330
      write(unitno,9010) rmmodes
 
331
c
 
332
      if (plot =='normal') then
 
333
         write(unitno,9011)
 
334
         do i=1,rmmodes             
 
335
            write(unitno,9012) i, rminfo(i,1), rminfo(i,4)
 
336
         end do
 
337
            write(unitno,9006)
 
338
      else if (plot =='resonance') then
 
339
         write(unitno,9013)
 
340
         do i=1,rmmodes
 
341
            write(unitno,9014) i, rminfo(i,1), rminfo(i,4), 
 
342
     &                        rminfo(i,2), rminfo(i,3)
 
343
         end do
 
344
            write(unitno,9006)
 
345
      end if ! normal
 
346
c
 
347
c ----------convert laser from au. -> nm -> 1/cm ----------------------------
 
348
c
 
349
        lambda = nm2icm/(au2nm/laser(1))
 
350
c
 
351
c ------------Calculate differential cross section-----------------------
 
352
c
 
353
      do i=start,rmmodes
 
354
         crs_real(i) = zero
 
355
         crs_imag(i) = zero
 
356
         crs_tot(i)  = zero
 
357
         tempfc=1.0D+06/(one-exp(-exparg*rminfo(i,1)/temp))
 
358
         frq4th = (lambda-rminfo(i,1))**4
 
359
         crs_tot(i) = tempfc*frq4th*conver*rminfo(i,4)/
 
360
     *                     (forty_five*rminfo(i,1))
 
361
         if (plot == 'resonance') then
 
362
            crs_real(i) = tempfc*frq4th*conver*rminfo(i,2)/
 
363
     *                     (forty_five*rminfo(i,1))
 
364
            crs_imag(i) = tempfc*frq4th*conver*rminfo(i,3)/
 
365
     *                     (forty_five*rminfo(i,1))
 
366
         end if
 
367
      enddo ! rmmodes
 
368
c
 
369
         frqmin = rminfo(1,1) - 60.0D+00  ! determine plot minimum
 
370
         frqmax = rminfo(rmmodes,1) + 60.0D+00 ! determine plot maximum
 
371
         dfrq   = ((frqmax-frqmin)) / (numpts-1) !energy steps determined by numpts
 
372
c
 
373
c ------calculate lineshape factors for the plot-----------------------
 
374
c
 
375
      if (line == 'lorentzian') then
 
376
         afac = width/(two*pi)
 
377
         bfac = width/two
 
378
      end if
 
379
c      
 
380
      if (line == 'gaussian') then
 
381
         width = width/widthfac
 
382
         afac = one/(sqrt(two*pi)*width)
 
383
         bfac = one/(two*width**two)
 
384
      end if
 
385
c
 
386
c ----------- make the plot---------------------------------------
 
387
c
 
388
      do i=1,numpts ! loop over number of points
 
389
         conv_tot = zero
 
390
         conv_real = zero
 
391
         conv_imag = zero
 
392
         frq = frqmin + i*dfrq
 
393
         do j=start,rmmodes ! loop over number of raman active modes
 
394
            if (line == 'lorentzian') then  
 
395
             factor = scale*afac/((frq-rminfo(j,1))**2 + bfac**2)
 
396
            else if (line == 'gaussian') then  
 
397
             factor = scale*afac*exp(-bfac* (frq-rminfo(j,1))**2)
 
398
            endif
 
399
            conv_tot = conv_tot + factor*crs_tot(j)
 
400
            if (plot == 'resonance') then
 
401
               conv_real = conv_real + factor*crs_real(j)
 
402
               conv_imag = conv_imag + factor*crs_imag(j)
 
403
            end if
 
404
         end do  ! end loop over rmmodes 
 
405
c
 
406
         if (plot == 'normal' ) then
 
407
            write (unitno,9017) frq, conv_tot
 
408
         else if (plot =='resonance') then
 
409
            write(unitno,9018) frq, conv_tot, conv_real, conv_imag
 
410
         endif  ! plot
 
411
      end do ! end loop over numpts
 
412
c
 
413
      close(unitno) 
 
414
c
 
415
      write(luout,22) filename(1:inp_strlen(filename))
 
416
 22   format(/' Raman scattering data written to ',a/)
 
417
      call util_flush(luout)
 
418
c
 
419
c --------------------------------------------------------------------
 
420
 9002 FORMAT(/1X,60(1H-)/5X,A16,'Raman Scattering Plot'/1X,60(1H-))
 
421
 9003 FORMAT(1X,'Excitation wavelength',F8.2,1X,'nm')
 
422
 9004 FORMAT(1X,'Convolution lineshape is ',A12)
 
423
 9005 FORMAT(1X,'FWHM',F8.2,1X,' 1/cm'/1X,60(1H-))
 
424
 9006 FORMAT(1X,60(1H-))
 
425
 9010 FORMAT(1X,60(1H-)/1X,'# of frequencies : ',I5)
 
426
 9011 FORMAT(1X,60(1H-)/1X,'#  freq [ 1/cm]  S [Ang**4/amu]'/
 
427
     *        1X,60(1H-))
 
428
 9012 FORMAT(1X,I3,F9.2,F14.4)
 
429
 9013 FORMAT(1X,60(1H-)/'#  freq [ 1/cm]  S-tot [a.u.]',  
 
430
     * 'S-real [a.u.]  S-imag [a.u.]'/1X,60(1H-))
 
431
 9014 FORMAT(1X,I3,F9.2,3F14.4)
 
432
 9015 FORMAT(1X,I4,F10.2,2E12.4)
 
433
 9016 FORMAT(1X,F10.2,3E12.4)
 
434
 9017 FORMAT(1X,F8.2,G14.4)
 
435
 9018 FORMAT(1X,F8.2,3G14.4)
 
436
c
 
437
      return
 
438
      end  
 
439
c ----------------end raman_scattering---------------------------------------
 
440
c
 
441
      subroutine raman_save(rtdb,ii,junk3,junk4)
 
442
c
 
443
      implicit none
 
444
#include "errquit.fh"
 
445
#include "rtdb.fh"
 
446
#include "mafdecls.fh"
 
447
#include "stdio.fh"
 
448
#include "geom.fh"
 
449
#include "nwc_const.fh" 
 
450
#include "inp.fh"
 
451
#include "global.fh"
 
452
c
 
453
      integer i,j,n,ii ! counting indexes
 
454
      integer rtdb     ! [input] rtdb handle
 
455
      double precision alfare(3,3),alfaim(3,3) ! [aoresponse] real and imaginary polarizanbility
 
456
      double precision junk3(6,3),junk4(6,3) ! temp. arrays for FD of polarizanbility
 
457
      logical debug
 
458
      logical status
 
459
c
 
460
      debug      = (.false. .and. ga_nodeid().eq.0) ! produce debug output
 
461
      status = rtdb_parallel(.true.)
 
462
c  -------------get raman ploarizability from rtdb ----------------------------
 
463
        if (.not. rtdb_get(rtdb, 'raman:alfare ', mt_dbl, 9, alfare))
 
464
     &   call errquit('raman_save:failed to get alfare', 1, RTDB_ERR)
 
465
        if (.not. rtdb_get(rtdb, 'raman:alfaim ', mt_dbl, 9, alfaim)) 
 
466
     &   call errquit('raman_save:failed to get alfaim', 2, RTDB_ERR)
 
467
c
 
468
      n=(ii-1)*3
 
469
      DO i=1,3
 
470
         n=n+1
 
471
        DO j=1,3
 
472
          junk3(n,j)= alfare(i,j) ! add AORESPONSE real data to array for finite_differencing
 
473
          junk4(n,j)= alfaim(i,j) ! add AORESPONSE imaginary data to array for finite_differencing
 
474
        enddo ! j loop
 
475
       enddo ! i loop
 
476
c  -------------delete raman ploarizability from rtdb ----------------------------
 
477
      IF (debug) THEN
 
478
          write(luout,*) 'junk3'
 
479
          call output(junk3,1,6,1,3,6,3,1)
 
480
          write(luout,*) 'junk4'
 
481
          call output(junk4,1,6,1,3,6,3,1)
 
482
          write(luout,*) 'alfare'
 
483
          call output(alfare,1,3,1,3,3,3,1)
 
484
          write(luout,*) 'alfaim'
 
485
          call output(alfaim,1,3,1,3,3,3,1)
 
486
      ENDIF !  for DEBUG
 
487
c
 
488
      return
 
489
      end
 
490
c ----------------end raman_save------------------------------------------------
 
491
c $Id: raman.F 21176 2011-10-10 06:35:49Z d3y133 $