From 79adcb335fa45b2abaa39d5daca38f2d56b31c6b Mon Sep 17 00:00:00 2001 From: edoapra Date: Tue, 8 Aug 2023 17:57:51 -0700 Subject: [PATCH] new geom keyword store_symrot to print gradients with the original geometry before autosym --- src/geom/geom_input.F | 23 ++++++++++++++++++++++- src/gradients/grad_force.F | 27 +++++++++++++++++++++++++-- 2 files changed, 47 insertions(+), 3 deletions(-) diff --git a/src/geom/geom_input.F b/src/geom/geom_input.F index edff6b9a0c..05e8b79bab 100644 --- a/src/geom/geom_input.F +++ b/src/geom/geom_input.F @@ -32,6 +32,7 @@ c logical ozmatrix logical oautoz logical oautosym + logical store_symrot logical ocenter logical found_cart logical ofinite, oatomfinite @@ -92,8 +93,10 @@ c logical found_ang2au_conv logical found_au2ang_conv double precision ang2au, au2ang + double precision rot(3,3),rot_inv(3,3) integer istart,iend,irate integer atn, isys, fcent + integer l_xyzold,k_xyzold character*16 format, nucmodel logical geom_zmtmak external geom_zmtmak @@ -122,7 +125,7 @@ c c c geometry [] [units ] [print] [noprint] [bqbq] \ c [angstrom_to_au ] [au_to_angstrom ] \ -c [(no)autoz] [nucleus ] +c [(no)autoz] [nucleus ] [store_symrot] c c subdirectives c @@ -161,6 +164,7 @@ c ecce_periodic = .false. sysset = .false. symset = .false. + store_symrot = .false. c call dfill(3*max_center, 0.0d0, velocities, 1) c @@ -257,6 +261,8 @@ c ------ Store ofinite to be used by HFine finite calc ---FA-03-21-11-- START $ call errquit('geo ofinite: rtdb_put failed',555, RTDB_ERR) c ------ Store ofinite to be used by HFine finite calc ---FA-03-21-11-- END end if + else if (inp_compare(.false.,'store_symrot',field)) then + store_symrot = .true. else if (name .ne. ' ') call errquit $ ('geom_input: geometry [] [units ]', 0, @@ -654,10 +660,25 @@ c write(luout,*) endif c + if(store_symrot) then + if(.not. ma_push_get(mt_dbl,3*ncenter,'xyzold', + L l_xyzold,k_xyzold)) call + $ errquit('geom_input: ma error', ncenter, MA_ERR) + call + D dcopy(ncenter*3, coords, 1, dbl_mb(k_xyzold), 1) + endif call geom_getsym0(rtdb,geom,oautosym, , ncenter,isys, , name,tags, , threquiv,scale,coords,charge,velocities,atomct) + if(store_symrot) then + call geom_findrot(ncenter,dbl_mb(k_xyzold),coords,rot) + call geom_m3inv(rot,rot_inv) + if (.not. rtdb_put(rtdb, 'geom:symrot', MT_DBL, 9, rot_inv)) + E call errquit("geom_input: symrot",0,0) + if (.not. ma_pop_stack(l_xyzold)) + & call errquit('geom_input:cant pop stack',0,0) + endif endif ! End of else clause for if(oadjust) c diff --git a/src/gradients/grad_force.F b/src/gradients/grad_force.F index 7fc2c5f245..26f7b986a1 100644 --- a/src/gradients/grad_force.F +++ b/src/gradients/grad_force.F @@ -168,14 +168,17 @@ c integer ifocc logical oskel, omp2, odft, ocdfit, status,frac_occ double precision xfac(numfunc), jfac, kfac + double precision rot(3,3),xrot,yrot,zrot c logical has_frac_occ external has_frac_occ c + double precision fol(3),f1,f2,f3 logical lcdft ! Qin integer ipop logical out1 logical zerof + logical have_symrot c cc AJL/Begin/SPIN ECPs integer ecp_channels @@ -233,7 +236,9 @@ c do_cosmo_smd = .false. endif endif - +c sym rot matrix + have_symrot=rtdb_get(rtdb, 'geom:symrot', MT_DBL, 9, rot) + C get SCF MO vectors for density odft = .false. @@ -1234,8 +1239,26 @@ C default: print the total forces do 30, i=1, nat if (.not. geom_cent_get(geom, i, tag, crd, q)) call errquit $ ('gradients: geometry corrupt?',0, GEOM_ERR) + if(have_symrot) then + xrot = crd(1)*rot(1,1) + crd(2)*rot(2,1) + crd(3)*rot(3,1) + yrot = crd(1)*rot(1,2) + crd(2)*rot(2,2) + crd(3)*rot(3,2) + zrot = crd(1)*rot(1,3) + crd(2)*rot(2,3) + crd(3)*rot(3,3) + crd(1)=xrot + crd(2)=yrot + crd(3)=zrot + fol(1)=dbl_mb(k_force+3*(i-1)) + fol(2)=dbl_mb(k_force+3*(i-1)+1) + fol(3)=dbl_mb(k_force+3*(i-1)+2) + f1 = fol(1)*rot(1,1) + fol(2)*rot(2,1) + fol(3)*rot(3,1) + f2 = fol(1)*rot(1,2) + fol(2)*rot(2,2) + fol(3)*rot(3,2) + f3 = fol(1)*rot(1,3) + fol(2)*rot(2,3) + fol(3)*rot(3,3) + else + f1=dbl_mb(k_force+3*(i-1)) + f2=dbl_mb(k_force+3*(i-1)+1) + f3=dbl_mb(k_force+3*(i-1)+2) + endif write(luout,2000) i, tag,(crd(j),j=1,3), - $ (dbl_mb(k_force+3*(i-1)+j),j=0,2) + $ f1,f2,f3 30 continue write(luout,*) 1000 format(/,/,25X,A,' ENERGY GRADIENTS',/,/,4X,'atom',15X,