1 subroutine qmd_forces(irtdb,isag,isat,xs,fs,msa,nsa,energy,lesp) 2c 3c $Id$ 4c 5 implicit none 6#include "errquit.fh" 7c 8#include "rtdb.fh" 9#include "geom.fh" 10#include "mafdecls.fh" 11#include "global.fh" 12#include "qmd_common.fh" 13#include "nwmd_constants.fh" 14c 15 external task_gradient,esp 16 logical task_gradient,esp 17c 18c data lgeom/0/ 19c 20 integer irtdb 21 integer msa,nsa 22 integer isag(msa),isat(msa) 23 real*8 xs(msa,3),fs(msa,3),energy 24 logical lesp 25c 26 integer i,isa 27 integer i_qmdq,l_qmdq,i_qmdt,l_qmdt,i_qmdb,l_qmdb 28 real*8 c(3),q 29 character*16 tag 30c 31 if(lgeom.eq.0) then 32 if(.not.geom_create(igeom,'geometry')) 33 + call errquit('qmd: failed to create geometry',0, GEOM_ERR) 34 if(.not.ma_push_get(mt_dbl,nsa,'qmdq',l_qmdq,i_qmdq)) 35 + call errquit('Failed to allocate qmdq',0, MA_ERR) 36 if(.not.ma_push_get(mt_byte,16*nsa,'qmdt',l_qmdt,i_qmdt)) 37 + call errquit('Failed to allocate qmdt',0, MA_ERR) 38 call qmd_geom_init(nsa,isag,isat,dbl_mb(i_qmdq),byte_mb(i_qmdt)) 39 if(ga_nodeid().eq.0) then 40 if(.not.geom_cart_set(igeom,nsa,byte_mb(i_qmdt),xs, 41 + dbl_mb(i_qmdq))) call errquit('Failed to initialize geom',0, 42 & GEOM_ERR) 43 endif 44 if(.not.geom_rtdb_store(irtdb,igeom,'geometry')) 45 + call errquit('Failed to store geom to rtdb',0, RTDB_ERR) 46 if(.not.ma_pop_stack(l_qmdt)) 47 + call errquit('Failed to deallocate qmdt',0, MA_ERR) 48 if(.not.ma_pop_stack(l_qmdq)) 49 + call errquit('Failed to deallocate qmdq',0, MA_ERR) 50 lgeom=1 51 endif 52c 53 if(.not.geom_rtdb_load(irtdb,igeom,'geometry')) 54 + call errquit('Failed to load geometry from rtdb',0, RTDB_ERR) 55c 56 do 1 i=1,nsa 57 isa=isag(i) 58 if(.not.geom_cent_get(igeom,isa,tag,c,q)) 59 + call errquit('geom_cent_get failed',isa, GEOM_ERR) 60 c(1)=xs(isa,1)*cnm2au 61 c(2)=xs(isa,2)*cnm2au 62 c(3)=xs(isa,3)*cnm2au 63 if(.not.geom_cent_set(igeom,i,tag,c,q)) 64 + call errquit('geom_cent_set failed',isa, GEOM_ERR) 65 1 continue 66 if(.not.geom_rtdb_store(irtdb,igeom,'geometry')) 67 + call errquit('geom_rtdb_store failed',0, RTDB_ERR) 68c 69 call ga_sync() 70 if(jdebug.gt.0) call ma_summarize_allocated_blocks() 71 call ga_sync() 72 if(.not.task_gradient(irtdb)) 73 + call errquit('qmd: Task_gradient failed',0, CALC_ERR) 74 call ga_sync() 75 if(lesp) then 76 if(.not.esp(irtdb)) call errquit('qmd esp error',0, RTDB_ERR) 77 endif 78c 79 if(.not.ma_push_get(mt_dbl,3*msa,'qmdb',l_qmdb,i_qmdb)) 80 + call errquit('Failed to allocate qmdq',0, MA_ERR) 81 if(.not.ma_verify_allocator_stuff()) 82 + call errquit('Error qmd_forces 1',0, MA_ERR) 83 call qmd_f(irtdb,isag,dbl_mb(i_qmdb),fs,msa,nsa,energy) 84 if(.not.ma_verify_allocator_stuff()) 85 + call errquit('Error qmd_forces 2',0, MA_ERR) 86 if(.not.ma_pop_stack(l_qmdb)) 87 + call errquit('Failed to deallocate qmdb',0, MA_ERR) 88c 89 return 90 end 91 subroutine qmd_f(irtdb,isag,buf,fs,msa,nsa,energy) 92c 93 implicit none 94#include "errquit.fh" 95c 96#include "rtdb.fh" 97#include "mafdecls.fh" 98#include "qmd_common.fh" 99#include "nwmd_constants.fh" 100c 101 integer irtdb,msa,nsa 102 integer isag(msa) 103 real*8 buf(3,msa),fs(msa,3),energy 104c 105 integer i,isa 106c 107 if(.not.rtdb_get(irtdb,'task:gradient',mt_dbl,3*msa,buf)) 108 + call errquit('qmd_f: rtdb_get gradient failed',0, RTDB_ERR) 109 if(.not.rtdb_get(irtdb,'task:energy',mt_dbl,1,energy)) 110 + call errquit('qmd_f: rtdb_get energy failed',0, RTDB_ERR) 111c 112 do 2 i=1,nsa 113 isa=isag(i) 114 fs(isa,1)=(-buf(1,i))*(cau2kj/cau2nm) 115 fs(isa,2)=(-buf(2,i))*(cau2kj/cau2nm) 116 fs(isa,3)=(-buf(3,i))*(cau2kj/cau2nm) 117 2 continue 118c 119 energy=energy*cau2kj 120c 121 return 122 end 123 subroutine qmd_geom_init(nsa,isag,isat,q,tag) 124c 125 implicit none 126c 127 integer cf_element 128 external cf_element 129c 130 integer nsa 131 integer isag(nsa),isat(nsa) 132 real*8 q(nsa) 133 character*16 tag(nsa) 134c 135 integer i,ia 136c 137 do 1 i=1,nsa 138 ia=cf_element(isat(i)) 139 call cf_num2tag(ia,tag(isag(i))) 140 q(i)=dble(ia) 141 1 continue 142c 143 return 144 end 145 subroutine qmd_wrttrj(lfntrj,mwm,nwm,mwa,nwa,msa,nsa, 146 + lxw,lvw,lxs,lvs,box,stime,pres,temp,tempw,temps,xw,vw,xs,vs) 147c 148 integer lfntrj,msa,nsa,mwm,nwm,mwa,nwa 149 logical lxw,lvw,lxs,lvs 150 real*8 stime,pres,temp,tempw,temps,box(3) 151 real*8 xw(mwm,3,mwa),vw(mwm,3,mwa) 152 real*8 xs(msa,3),vs(msa,3) 153c 154 character*10 rdate,rtime 155c 156#include "global.fh" 157c 158 if(ga_nodeid().eq.0) then 159c 160 call swatch(rdate,rtime) 161c 162 write(lfntrj,1000) 163 1000 format('frame') 164 write(lfntrj,1001) stime,temp,pres,rdate,rtime 165 1001 format(2f12.6,1pe12.5,1x,2a10) 166 write(lfntrj,1002) box(1),0.0d0,0.0d0 167 write(lfntrj,1002) 0.0d0,box(2),0.0d0 168 write(lfntrj,1002) 0.0d0,0.0d0,box(3) 169 1002 format(3f12.6) 170 write(lfntrj,1003) lxw,lvw,lxs,lvs,nwm,nwa,nsa 171 1003 format(4l1,3i10) 172c 173 if((lxw.or.lvw).and.nwm.gt.0) then 174 if(lvw) then 175 do 1 i=1,nwm 176 write(lfntrj,1004) ((xw(i,j,k),j=1,3),(vw(i,j,k),j=1,3),k=1,nwa) 177 1004 format(6f8.3) 178 1 continue 179 else 180 do 2 i=1,nwm 181 write(lfntrj,1005) ((xw(i,j,k),j=1,3),k=1,nwa) 182 1005 format(3f8.3) 183 2 continue 184 endif 185 endif 186c 187 if((lxs.or.lvs).and.nsa.gt.0) then 188 if(lvs) then 189 do 3 i=1,nsa 190 write(lfntrj,1006) (xs(i,j),j=1,3),(vs(i,j),j=1,3) 191 1006 format(6f8.3) 192 3 continue 193 else 194 do 4 i=1,nsa 195 write(lfntrj,1007) (xs(i,j),j=1,3) 196 1007 format(3f8.3) 197 4 continue 198 endif 199 endif 200c 201 endif 202c 203 return 204 end 205