1 subroutine vib_thermo(rtdb,nat,nat3, 2 & eigvals,masses,coords,scr) 3 implicit none 4#include "errquit.fh" 5* 6* $Id$ 7* 8* routine to compute the thermodynamic properties under the 9* rigid rotor harmonic approximation: 10* 11* See: 12* 1. "Statistical Mechanics" by Donald A. McQuarrie, Chapters 5,6,8, 13* Publisher Harper and Row, (c) 1976 14* 2. 15* 16*::functions:: 17#include "rtdb.fh" 18#include "mafdecls.fh" 19#include "util.fh" 20#include "stdio.fh" 21#include "inp.fh" 22*::passed:: 23 integer rtdb ! [input] rtdb handle 24 integer nat ! [input] number of atoms 25 integer nat3 ! [input] 3*nat 26 double precision eigvals(nat3) ! [input] eigenvalues (in cm**-1) 27 double precision masses(nat) ! [input] mass of atoms (in amu) 28 double precision coords(3,nat) ! [input] atomic coordinates (in a.u.) 29 double precision scr(*) 30*::local:: 31 integer i, j, count, error 32 integer iat 33 integer len_theory 34 integer nopen, multiplicity 35 logical linear 36 double precision com(3) ! center of mass coordinates 37 double precision Inertia(3,3), Inertia2(3,3), rx, ry, rz 38 double precision VdI(3,3) 39 double precision Ivals(3), Ivecs(3,3), Itri(3*4/2) 40 double precision total_mass, xmass 41 double precision H_trans , G_trans , S_trans , Cp_trans , Cv_trans 42 double precision H_rot , G_rot , S_rot , Cp_rot , Cv_rot 43 double precision H_vib , G_vib , S_vib , Cp_vib , Cv_vib 44 double precision H_elec , G_elec , S_elec , Cp_elec , Cv_elec 45 double precision H_total , G_total , S_total , Cp_total , Cv_total 46 double precision Temp 47 double precision R_gas 48 double precision Pi 49 double precision Sigma 50 character*32 theory 51 character*50 my_string 52* 53 double precision thresh 54 parameter (thresh=1.0d-3) 55c::-statement function 56 logical is_it_close_to 57 double precision value,test 58 intrinsic abs 59*--- is value close to test? 60 is_it_close_to(value,test) = (abs(value-test).lt.thresh) 61* 62* -- until it is working 63 if (nat.gt.0) return 64* 65 Pi = 2.0d00*acos(0.0d00) 66 H_trans = 0.0d00 67 H_rot = 0.0d00 68 H_vib = 0.0d00 69 H_elec = 0.0d00 70 H_total = 0.0d00 71 G_trans = 0.0d00 72 G_rot = 0.0d00 73 G_vib = 0.0d00 74 G_elec = 0.0d00 75 G_total = 0.0d00 76 S_trans = 0.0d00 77 S_rot = 0.0d00 78 S_vib = 0.0d00 79 S_elec = 0.0d00 80 S_total = 0.0d00 81 Cp_trans = 0.0d00 82 Cp_rot = 0.0d00 83 Cp_vib = 0.0d00 84 Cp_elec = 0.0d00 85 Cp_total = 0.0d00 86 Cv_trans = 0.0d00 87 Cv_rot = 0.0d00 88 Cv_vib = 0.0d00 89 Cv_elec = 0.0d00 90 Cv_total = 0.0d00 91*-- electronic component 92 my_string = ' ' 93 theory = ' ' 94 if (.not.rtdb_cget(rtdb, 'task:theory', 1, theory)) call errquit 95 & ('vib_thermo: could not read task:theory',911, RTDB_ERR) 96 len_theory = inp_strlen(theory) 97 my_string = theory(1:len_theory) // ':nopen' 98 nopen = 0 99 if (.not.rtdb_get(rtdb,my_string,mt_int,1,nopen)) then 100 write(luout,*)' theory is :',theory(1:len_theory) 101 call errquit('vib_thermo: could not read --theory--:nopen',911, 102 & RTDB_ERR) 103 endif 104 multiplicity = nopen + 1 105 S_elec = R_gas*log(dble(multiplicity)) 106 H_elec = 0.0d00 107 Cp_elec = 0.0d00 108 Cv_elec = 0.0d00 109 G_elec = H_elec - Temp * S_elec 110*-- rotational component 111*- compute total mass and center of mass 112 call dfill(3,0.0d00,com,1) 113 total_mass = 0.0d00 114 do iat = 1,nat 115 xmass = masses(iat) 116 total_mass = total_mass + xmass 117 com(1) = com(1) + xmass*coords(1,iat) 118 com(2) = com(2) + xmass*coords(2,iat) 119 com(3) = com(3) + xmass*coords(3,iat) 120 enddo 121 call dscal(3,(1.0d00/total_mass),com,1) 122* 123 write(luout,*)' Coordinates ' 124 call output(coords,1,3,1,3,3,3,1) 125 write(luout,*)' Center of Mass :',com 126* 127 call dfill((3*3),0.0d00,Inertia,1) 128 do iat = 1,nat 129 xmass = masses(iat) 130 rx = coords(1,iat) - com(1) 131 ry = coords(2,iat) - com(2) 132 rz = coords(3,iat) - com(3) 133 Inertia(1,1) = Inertia(1,1) + 134 & xmass*(ry*ry + rz*rz) 135 Inertia(2,1) = Inertia(2,1) - 136 & xmass*ry*rx 137 Inertia(3,1) = Inertia(3,1) - 138 & xmass*rz*rx 139 Inertia(2,2) = Inertia(2,2) + 140 & xmass*(rx*rx + rz*rz) 141 Inertia(3,2) = Inertia(3,2) - 142 & xmass*rz*ry 143 Inertia(3,3) = Inertia(3,3) + 144 & xmass*(rx*rx + ry*ry) 145 enddo 146 Inertia(1,2) = Inertia(2,1) 147 Inertia(1,3) = Inertia(3,1) 148 Inertia(2,3) = Inertia(3,2) 149 write(luout,*)' Raw Inertial Matrix' 150 call output(Inertia,1,3,1,3,3,3,1) 151* set up triangular matrix 152 count = 0 153 do i = 1,3 154 do j = 1,i 155 count = count + 1 156 Itri(count) = Inertia(i,j) 157 enddo 158 enddo 159* diagonilze Inertia Matrix (parinoia) 160 call dfill(3,0.0d00,Ivals,1) 161 call dfill((3*3),0.0d00,Ivecs,1) 162 call vib_sjacobi(3,3,Itri,Ivals,Ivecs,error) 163 if (error.ne.0) call errquit('vib_thermo: vib_sjacobi failed',911, 164 & UNKNOWN_ERR) 165* 166 write(luout,*) 'eigenvalues of Inertia Matrix' 167 call output(Ivals,1,3,1,1,3,1,1) 168 write(luout,*) 'eigenvectors of Inertia Matrix' 169 call output(Ivecs,1,3,1,3,3,3,1) 170* form "diagonal" matrix 171 call dfill((3*3),0.0d00,Inertia2,1) 172 call dgemm('t','n',3,3,3,1.0d00,Ivecs,3,Inertia,3,0.0d00,VdI,3) 173 call dgemm('n','n',3,3,3,1.0d00,VdI,3,Ivecs,3,0.0d00,Inertia2,3) 174 write(luout,*)' Diagonalized Inertial Matrix' 175 call output(Inertia2,1,3,1,3,3,3,1) 176* 177 linear = is_it_close_to(Ivals(1),0.0d00) .and. 178 & is_it_close_to((abs(Ivals(2)-Ivals(3))),0.0d00) 179* 180 sigma = 1 181 if (linear) then 182 S_rot = R_gas*(1 + log(Ivals(2)/Sigma*1.0d00)) 183 if (rtdb_put(rtdb,'vib:linear',mt_log,1,.true.)) 184 & call errquit('vib_thermo: failed to set linear',555, RTDB_ERR) 185 else 186 endif 187* 188 end 189