1C $Id$ 2 Subroutine nwc_gen (x,v,dx,fmat,hess,xs,natom,n3,n3tm,ndbar,op) 3C 4C Run NWChem task to get potential and first and second derivatives 5C for geometry x. 6C 7 Implicit None 8#include "errquit.fh" 9#include "nwc_drdyP.fh" 10#include "geom.fh" 11#include "mafdecls.fh" 12#include "printlevels.fh" 13#include "rtdb.fh" 14#include "stdio.fh" 15 Integer natom, n3, n3tm, ndbar, op 16 Double Precision v 17 Double Precision x(n3),dx(n3),fmat(n3tm,n3), 18 & hess(ndbar),xs(n3), velocities(3,n3) 19* 20 Integer ndima 21 Parameter (ndima=3) 22 Double Precision amat(ndima,ndima), 23 & rr(ndima,ndima),rrs(ndima,ndima), 24 * rri(ndima,ndima) 25* 26 Integer geom 27 Integer i,j,jj 28 double precision threquiv 29 character*16 groupname 30 character*255 dummy_file 31 Integer ncenter_B4_autosym, ncenter, nata 32 Integer current_print_level,nops 33 logical oautosym 34 logical status_rtdb, status_ignore 35* logical copy_sym, copy_c1, oautosym 36* 37 38 call util_print_get_level(current_print_level) 39 call util_print_set_level(print_none) 40* 41 if (op.eq.DRDY_CODE_HESSIAN) then 42* 43* ... delete finite difference files that may exist. 44 call util_file_name('hess', .false., .false.,dummy_file) 45 call util_file_unlink(dummy_file) 46 call util_file_name('fd_ddipole',.false.,.false.,dummy_file) 47 call util_file_unlink(dummy_file) 48 endif 49* 50 if (.not.geom_create(geom,'geometry')) 51 & call errquit('nwc_gen: geom_create failed',911, 52 & GEOM_ERR) 53 if (.not.geom_set_user_units(geom,'a.u.')) 54 & call errquit('nwc_gen: geom_set_user_units failed',911, 55 & GEOM_ERR) 56 if (.not.geom_cart_set(geom,natom,atomic_labels,x,atomic_charge)) 57 & call errquit('nwc_gen: geom_cart_set failed',911, 58 & GEOM_ERR) 59 if (.not.geom_masses_set(geom,natom,nwcmass)) 60 & call errquit('nwc_gen:geom_masses_set failed',911, 61 & GEOM_ERR) 62 if (.not.geom_atomct_set(geom,natom,atomic_atomct)) 63 & call errquit('nwc_gen:geom_atomct_set failed',911, 64 & GEOM_ERR) 65* if (.not.geom_rtdb_store(my_rtdb,geom,'geometry')) 66* & call errquit('nwc_gen: geom_rtdb_store failed',911, 67* & RTDB_ERR) 68c 69c Now figure out if we are going to use autosym and if so, get the 70c value for the threshold and take care of the symmetry 71c 72 status_rtdb = rtdb_parallel(.false.) 73 74 if (.not. rtdb_get(my_rtdb,'drdy:autosym',mt_log,1,oautosym)) 75 & call errquit('nwc_gen: problem getting autosym from rtdb',555, 76 & RTDB_ERR) 77 if (oautosym) then 78 if (.not.rtdb_get(my_rtdb,'drdy:threquiv',mt_dbl,1,threquiv)) 79 & call errquit('nwc_gen:problem getting threquiv from rtdb',555, 80 & RTDB_ERR) 81 ncenter_B4_autosym = natom 82 ncenter = natom 83 call dcopy(n3,x,1,dx,1) ! temporary use of dx 84 call dcopy(ncenter,atomic_charge,1,copy_charge,1) 85 call dcopy(ncenter,atomic_atomct,1,copy_atomct,1) 86 do i = 1,ncenter 87 copy_labels(i) = atomic_labels(i) 88 enddo 89 call dcopy(3*n3,0.d0,0,velocities,1) 90 call geom_auto_sym(my_rtdb,geom,dx, 91 & copy_charge,copy_labels,copy_atomct,ncenter, 92 & threquiv,groupname,velocities) 93 if (op.ne.DRDY_CODE_SPENERGY) 94 & write(luout,10) groupname 95 if (groupname(1:2).ne."C1") then 96 if (geom_group_set(geom,groupname)) then 97 if (.not.geom_cart_set(geom,ncenter,copy_labels,dx, 98 & copy_charge)) 99 & call errquit('nwc_gen: geom_cart_set failed', 0, 100 & GEOM_ERR) 101 if (ncenter_B4_autosym .ne. ncenter) call errquit 102 & ('nwc_gen: autosym bug : number of atoms wrong', 103 & ncenter, INPUT_ERR) 104 else 105 write(luout,*) ' autosym detected unknown group ', 106 & groupname 107 call errquit('nwc_gen: autosym: invalid group',0, INPUT_ERR) 108 endif 109c 110c Apply system and symmetry info to the list of 111c unique centers build mapping tables set up coord lists 112c 113 nata = ncenter 114 call sym_nwc(geom,my_rtdb,nata,.false.,1.0d00,threquiv,nops, 115 % .false.) 116c 117c Check that if we used autosym that we ended up with the 118c same no. of atoms ... if we don't then autosym and nwchemsym 119c don't agree on the orientation of point group elements 120c 121 if (.not. geom_ncent(geom,ncenter)) call errquit 122 $ ('nwc_gen: geom_cent?',0, GEOM_ERR) 123 if (ncenter_B4_autosym .ne. ncenter) call errquit 124 $ ('nwc_gen: autosym bug : too many atoms',ncenter, 125 & INPUT_ERR) 126c 127c Force exact symetry on the coordinates 128c 129 call sym_geom_project(geom, threquiv) 130 endif 131 endif 132c 133c store as default geometry 134c 135 if (.not.geom_rtdb_store(my_rtdb,geom,'geometry')) 136 & call errquit('nwc_gen: geom_rtdb_store failed',911, 137 & RTDB_ERR) 138 if (.not.geom_destroy(geom)) 139 & call errquit('nwc_gen: geom_destroy failed',911, GEOM_ERR) 140 call util_print_set_level(current_print_level) 141* write(luout,*)' x after ',x 142c 143c Actually do the deed 144c 145 status_ignore = rtdb_parallel(status_rtdb) 146 status_rtdb = rtdb_parallel(.false.) 147 148 if ((op.ge.DRDY_CODE_SPENERGY).and. 149 & (op.le.DRDY_CODE_HESSIAN)) then 150 call drdy_synch(op,'nwc_gen') 151 else 152 write(luout,*) 'Trying to calculate an unknown op code' 153 call errquit('nwc_gen: Unknown op code',op, INPUT_ERR) 154 endif 155 status_ignore = rtdb_parallel(status_rtdb) 156 status_rtdb = rtdb_parallel(.false.) 157c 158c Get the appropriate values to hand back 159c 160 call drdy_nwc_get_energy(my_rtdb,v) 161 if (op.ge.DRDY_CODE_GRADIENT) then 162 call drdy_nwc_get_coords(my_rtdb,xs) 163 call drdy_nwc_get_gradient(my_rtdb,dx) 164 endif 165 if (op.eq.DRDY_CODE_HESSIAN) 166 & call drdy_nwc_get_hessian(my_rtdb,hess) 167* 168C 169C Get transformation matrix (rotation matrix amat) from the gaussian 170C standard orientation (xs) to the original orientation (x) 171 if (op.ge.DRDY_CODE_GRADIENT) then 172 call drdy_rotmat(x,xs,amat,rr,rrs,rri,natom,ndima) 173 call drdy_rotg(dx,amat,rr,natom,n3) 174 endif 175C 176 if (op.eq.DRDY_CODE_HESSIAN) then 177 jj = 0 178 do i = 1,n3 179 do j = 1,i 180 jj = jj + 1 181 fmat(i,j) = hess(jj) 182 fmat(j,i) = hess(jj) 183 enddo 184 enddo 185C Transform hessian matrix from gaussian standard orientation to the 186C original orientation 187 call drdy_rotf(fmat,amat,rr,natom,n3,n3tm) 188 endif 189* 190 status_ignore = rtdb_parallel(status_rtdb) 191 return 19210 Format (1x,'Symmetry group is ',a3) 1931000 Format (1x,79a1) 1941001 Format (' ') 1951002 Format (1x,2i2) 1961003 Format (1x,5a1,1p3e20.10) 197 End 198