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