From bf72729c8983af05c6552e6185dfa803b4e1dc86 Mon Sep 17 00:00:00 2001 From: Sean Fischer Date: Mon, 1 Feb 2016 18:48:44 +0000 Subject: [PATCH] fixed the size of the step that was taken for finite difference, fixed variables so calculations on linear molecules now work --- src/property/raman.F | 31 +++++++++++------- src/property/raman_input.F | 8 ++--- src/property/task_raman.F | 64 +++++++++++++++++++++++++++----------- 3 files changed, 69 insertions(+), 34 deletions(-) diff --git a/src/property/raman.F b/src/property/raman.F index 4b1b863d69..0c63e8e8f1 100644 --- a/src/property/raman.F +++ b/src/property/raman.F @@ -29,8 +29,8 @@ c integer rtdb ! [input] rtdb handle integer natom ! [input] number of atoms integer nat3 ! [input] 3*number of atoms - integer first ! first mode to consider in aoresponse (default =7 ramana =1 hyperraman) - integer tmpmode ! set to fill rminfo from 1 ( not 7 for raman calc) + integer first ! first mode to consider in aoresponse (default =6 or 7 raman =1 hyperraman) + integer tmpmode ! set to fill rminfo from 1 ( not 6 or 7 for raman calc) integer rmmodes ! # of raman active modes double precision rminfo(rmmodes,4) ! data for raman spec @@ -40,6 +40,10 @@ c double precision eigenvals(nat3) ! [input] (mode) double precision ncoords(3,natom) ! [scratch] coords after step double precision steps(3,natom) ! [scratch] step generated by vector and scaled +c + double precision length_of_step, scale + double precision ddot + external ddot c parameter (bohr2ang=0.52917724924D+00) ! CONVERSION OF BOHR to ANGSTROMS c -------------determine sign of the step--------------------------------- @@ -57,13 +61,16 @@ c 1 format('FA-check-rminfo:(imode,iii,first,tmpmode,rminfo)=(', c & i4,',',i4,',',i4,',',i4,',',f15.8,')') c ======= FA-check rminfo(x,1) ======== END c -------------------------------------------------------------------- - ivec = 1 - do iatom = 1,natom - do ixyz = 1,3 - steps(ixyz,iatom)=sign*step_size*eigenvecs(ivec,imode) - ivec = ivec + 1 - enddo ! ixyz - enddo ! iatom + ivec = 1 + do iatom = 1,natom + do ixyz = 1,3 + steps(ixyz,iatom)=eigenvecs(ivec,imode) + ivec = ivec + 1 + enddo ! ixyz + enddo ! iatom + length_of_step = sqrt(ddot(nat3,steps,1,steps,1)) + scale = sign*step_size/length_of_step + call dscal(nat3,scale,steps,1) call daxpy(nat3,1.0d00,steps,1,ncoords,1) ! mult coords if (.not. geom_cart_coords_set(geom,ncoords)) @@ -85,7 +92,8 @@ c -------------------------------------------------------------------- & step_size,! in : step of finite differencing & rminfo, ! in : Raman data & eigenvecs,! in : normal modes eigenvectors (nat3,nat3) - & mass) ! in : mass + & mass, ! in : mass + & first0) ! in : first nonzero mode (6 or 7) c c Authors: Jonathan Mullin, Northwestern University (ver 1: Jan. 2011) c Fredy W. Aquino, Northwestern University (ver 2: Oct. 2012) @@ -108,6 +116,7 @@ c integer imode ! mode # integer natom ! [input] number of atoms integer nat3 ! [input] 3*number of atoms + integer first0 ! [input] first nonzero mode (6 or 7) c double precision rminfo(rmmodes,4) ! raman data double precision step_size,stepsize ! [input] step of finite differencing @@ -134,7 +143,7 @@ c -------------------------------------------------------------------- call dfill(3*natom,0.0D+00,tmode,1) ! c zero stepsize = zero - m = imode - 6 + m = imode - first0 + 1 j=1 i=1 ar2 = zero ! alpha real diff --git a/src/property/raman_input.F b/src/property/raman_input.F index bcd035a61f..1a9972d2b0 100644 --- a/src/property/raman_input.F +++ b/src/property/raman_input.F @@ -47,6 +47,7 @@ c ---------- FA-definitions -------- END c c set some defaults c + field=' ' plot = 'normal' ! normal or resonance line = 'lorentzian' ! lorentzian (l) or gaussian (g) lineshape width = 20.0D+00 ! full-width at half maximum (FWHM) in 1/cm @@ -54,7 +55,6 @@ c hyperraman = .false. ! flag to calculate hyperaman terms vroa = .false. ! flag to calculate vibrational raman spec rmmodes = 0 - first = 7 last = 10000 low = 0.0D+00 high = 100000.0D+00 @@ -132,9 +132,9 @@ c --- determine first normal mode to use --- else if(inp_compare(.false.,'first',test)) then if(.not. inp_i(first)) $ call errquit(pname//'missing value for first',911, INPUT_ERR) - if (.not. rtdb_put(rtdb,'raman:first',mt_int,1,first)) - $ call errquit(pname//'rtdb put failed',0, RTDB_ERR) -c --- determine first normal mode to use --- +c --- not setting default here, it will be set later after +c frequency calculation has been done so we know if we have +c a linear molecule or not else if(inp_compare(.false.,'last',test)) then if(.not. inp_i(last)) ! FA-06-16-12 bug-fixed (BEF: first AFT: last) $ call errquit(pname//'missing value for last',911, INPUT_ERR) diff --git a/src/property/task_raman.F b/src/property/task_raman.F index 15ebf284eb..2711d6933e 100644 --- a/src/property/task_raman.F +++ b/src/property/task_raman.F @@ -59,6 +59,7 @@ c --> Experimental (not published yet) integer j,pos,first0 ! FA-06-15-12 logical preraman ! FA-06-18-12 + logical linear character*32 pname @@ -107,6 +108,12 @@ c endif $ call errquit(pname//'rtdb_put freq_done',911, RTDB_ERR) endif c +c --------Figure out if molecule is linear------------ + +c if vib module doesn't list molecule as linear, assume it is not + if (.not. rtdb_get(rtdb,'vib:linear',mt_log,1,linear)) + $ linear=.false. +c c --------Create/load reference geometry to get the number of atoms------------ if (.not.geom_create(geom,'geometry')) call errquit @@ -116,7 +123,11 @@ c --------Create/load reference geometry to get the number of atoms------------ if (.not. geom_ncent(geom,nat)) & call errquit(pname//'geom_ncent failed?',3, GEOM_ERR) nc = nat*3 - rmmodes = nc-6 + if (linear) then + rmmodes = nc-5 + else + rmmodes = nc-6 + end if c if (ga_nodeid().eq.0) then c write(*,1) nat,nc,rmmodes @@ -146,8 +157,13 @@ c ----------Load parameters from rtdb or assign defaults----------------------- $ low = 0.0D+00 ! lowest wavenumber normal mode to use if (.not. rtdb_get(rtdb,'raman:high',mt_dbl,1,high)) $ high = 100000.0D+00 ! Highest wavenumber normal mode to use - if (.not. rtdb_get(rtdb,'raman:first',mt_int,1,first)) - $ first = 7 ! first normal mode to use + if (.not. rtdb_get(rtdb,'raman:first',mt_int,1,first)) then + if (linear) then + first = 6 ! first normal mode to use + else + first = 7 ! first normal mode to use + end if + end if if (.not. rtdb_get(rtdb,'raman:last',mt_int,1,last)) $ last = 10000 ! last normal mode to use if (.not. rtdb_get(rtdb,'raman:hyperraman',mt_log,1,hyperraman)) @@ -156,7 +172,11 @@ c ----------Load parameters from rtdb or assign defaults----------------------- $ vroa = .false. ! # flag to calculate vibrational if (.not. rtdb_get(rtdb,'raman:preraman',mt_log,1,preraman)) $ preraman = .false. ! # flag to do task_freq() and leave - first0=7 ! constant + if (linear) then + first0=6 ! constant + else + first0=7 ! constant + end if c ======== FA-debug =============== START c if (ga_nodeid().eq.0) then c write(*,2) plot,line,width,step_size,steps @@ -172,8 +192,13 @@ c ======== FA-debug =============== END rmmodes = nc c c --- in case we want overide the defaults for modes to include --- - if (.not. rtdb_get(rtdb,'raman:first',mt_int,1,first)) - $ first = 7 ! srtep size for displacement along modes + if (.not. rtdb_get(rtdb,'raman:first',mt_int,1,first)) then + if (linear) then + first = 6 ! srtep size for displacement along modes + else + first = 7 ! srtep size for displacement along modes + end if + end if endif c c ----------alocate space for freq and normal modes---------------------------- @@ -294,7 +319,7 @@ c c ------------enough setup really do the calculation------------------------ if (.not.preraman) then call task_raman_doit(rtdb,geom,nc,nat, - & first0, ! = 7 constant + & first0, ! = 6 or 7 & first,last,rmmodes, & steps,nfreq,plot,line,width, & step_size, @@ -336,7 +361,7 @@ c c c == perform raman calculation == subroutine task_raman_doit(rtdb,geom,nc,nat, - & first0, ! = 7 constant + & first0, ! = 6 or 7 & first,last,rmmodes, & steps,nfreq, & plot,line,width, @@ -495,7 +520,7 @@ c 0000000 Reading modes from input script 000000000 START & lbl_raman, ! in: raman label & begin, ! in: & last, ! in: - & first0, ! in: = 7 constant + & first0, ! in: = 6 or 7 & eigenvecs, ! in: hessian data (modes) & eigenvals, ! in: hessian data (frequencies) & mass, ! in: mass(i) i=1,nat @@ -519,7 +544,7 @@ c 0000000 Reading modes from input script 000000000 START & lbl_raman, ! in: raman label & mode_ini, ! in: & mode_end, ! in: - & first0, ! in: = 7 constant + & first0, ! in: = 6 or 7 & eigenvecs, ! in: hessian data (modes) & eigenvals, ! in: hessian data (frequencies) & mass, ! in: mass(i) i=1,nat @@ -541,7 +566,7 @@ c 0000000 Reading modes from input script 000000000 START & lbl_raman, ! in: raman label & begin, ! in: starting mode & last, ! in: ending mode - & first0, ! in: = 7 constant + & first0, ! in: = 6 or 7 & eigenvecs, ! in: hessian data (modes) & eigenvals, ! in: hessian data (frequencies) & mass, ! in: mass(i) i=1,nat @@ -596,7 +621,7 @@ c 00000000000000000000000000000000000000000000000000000000000000 & rmmodes, ! in: total nr. modes & rminfo, ! in: stores raman info & nc,nat, ! in: (nc,nat)=(nr coord,nr atoms) - & first0, ! in: = 7 constant + & first0, ! in: = 6 or 7 & eigenvecs, ! in: hessian data (modes) & eigenvals, ! in: hessian data (frequencies) & mass, ! in: mass(i) i=1,nat @@ -757,7 +782,8 @@ c ------------------------do finite defernced for mode ii---------------------- & step_size, & rminfo, & eigenvecs, - & mass) + & mass, + & first0) c ======== FA: Writing to file rminfo ========= START c if (ga_nodeid().eq.0) c & write(*,*) 'BEF raman_write() ...' @@ -783,7 +809,7 @@ c 00000000000000000000000000000000000000000000000000000000000000 & lbl_raman, ! in: raman label & begin, ! in: starting mode & last, ! in: ending mode - & first0, ! in: = 7 constant + & first0, ! in: = 6 or 7 & eigenvecs, ! in: hessian data (modes) & eigenvals, ! in: hessian data (frequencies) & mass, ! in: mass(i) i=1,nat @@ -890,7 +916,7 @@ c ======== construct filename in listing mode ====== END & rmmodes, ! in: total nr. modes & rminfo, ! in: stores raman info & nc,nat, ! in: (nc,nat)=(nr coord,nr atoms) - & first0, ! in: = 7 constant + & first0, ! in: = 6 or 7 & eigenvecs, ! in: hessian data (modes) & eigenvals, ! in: hessian data (frequencies) & mass, ! in: mass(i) i=1,nat @@ -915,7 +941,7 @@ c ======== construct filename in listing mode ====== END & lbl_raman, ! in: raman label & mode_ini, ! in: & mode_end, ! in: - & first0, ! in: = 7 constant + & first0, ! in: = 6 or 7 & eigenvecs, ! in: hessian data (modes) & eigenvals, ! in: hessian data (frequencies) & mass, ! in: mass(i) i=1,nat @@ -1036,7 +1062,7 @@ c ************************************************************ & rmmodes, ! in: total nr. modes & rminfo, ! in: stores raman info & nc,nat, ! in: (nc,nat)=(nr coord,nr atoms) - & first0, ! in: = 7 constant + & first0, ! in: = 6 or 7 & eigenvecs, ! in: hessian data (modes) & eigenvals, ! in: hessian data (frequencies) & mass, ! in: mass(i) i=1,nat @@ -1058,7 +1084,7 @@ c ************************************************************ & lbl_raman, ! in: raman label & begin, ! in: & last, ! in: - & first0, ! in: = 7 constant + & first0, ! in: = 6 or 7 & eigenvecs, ! in: hessian data (modes) & eigenvals, ! in: hessian data (frequencies) & mass, ! in: mass(i) i=1,nat @@ -1139,7 +1165,7 @@ c if (begin .gt. last) go to 200 & rmmodes, ! in: total nr. modes & rminfo, ! in: stores raman info & nc,nat, ! in: (nc,nat)=(nr coord,nr atoms) - & first0, ! in: = 7 constant + & first0, ! in: = 6 or 7 & eigenvecs, ! in: hessian data (modes) & eigenvals, ! in: hessian data (frequencies) & mass, ! in: mass(i) i=1,nat