1*
2* $Id$
3*
4
5*     ***************************
6*     *                         *
7*     *    seperate_molpsp      *
8*     *                         *
9*     ***************************
10      subroutine seperate_molpsp(rtdb)
11      implicit none
12#include "errquit.fh"
13      integer rtdb
14
15#include "bafdecls.fh"
16#include "btdb.fh"
17#include "geom.fh"
18
19
20
21      logical bqbq
22      common / nwpw_bqbq_block / bqbq
23
24      logical mmexist
25      common / ion_qmmm/ mmexist
26
27      integer nionall
28      common / ion_nionall_block / nionall
29
30*     *** local variables ***
31      integer     i,geom1,geom2,geom3,nion1,nion2,nion3
32      integer rt(2),tt(2),qt(2),mt(2)
33      double precision q,rxyz(3)
34      character*16     t
35      logical value
36
37*     **** external functions ****
38      logical          parseqmmm
39      external         parseqmmm
40
41
42*     *******************************************************************
43*     **** seperate ions and charges from molecular pseudopotentials ****
44*     *******************************************************************
45
46*     **** generate chargepspwgeometry from geometry ****
47*     **** generate qmmmgeometry from geometry ****
48      value = geom_create(geom1,'geometry')
49      value = value.and.geom_create(geom2,'chargepspwgeometry')
50      value = value.and.geom_create(geom3,'qmmmgeometry')
51      value = value.and.geom_rtdb_load(rtdb,geom1,'geometry')
52      value = value.and.geom_ncent(geom1,nion1)
53      if (.not. value) call errquit('opening geometry',0, GEOM_ERR)
54
55*     *** set nionall and bqbq ****
56      nionall = nion1
57      bqbq    = geom_include_bqbq(geom1)
58
59      value = BA_push_get(mt_dbl, (3*nion1), 'rt',rt(2),rt(1))
60      value = value.and.
61     >        BA_push_get(mt_dbl, (nion1),   'qt',qt(2),qt(1))
62      value = value.and.
63     >        BA_push_get(mt_dbl, (nion1),   'mt',mt(2),mt(1))
64      value = value.and.
65     >        BA_push_get(mt_byte,(16*nion1),'tt',tt(2),tt(1))
66      if (.not. value) call errquit('out of stack memory',0, MA_ERR)
67
68      nion2 = 0
69      nion3 = 0
70      do i=1,nion1
71         value = geom_cent_get(geom1,i,t,rxyz,q)
72         if (.not.parseqmmm(t)) then
73            nion2 = nion2 + 1
74         else
75            nion3 = nion3 + 1
76         end if
77      end do
78      value = value.and.
79     >        geom_cart_get(geom1,nion1,byte_mb(tt(1)),
80     >                                  dbl_mb(rt(1)),
81     >                                  dbl_mb(qt(1)))
82      value = value.and.
83     >        geom_cart_set(geom2,nion2,byte_mb(tt(1)),
84     >                                  dbl_mb(rt(1)),
85     >                                  dbl_mb(qt(1)))
86      value = value.and.geom_masses_get(geom1,nion1,dbl_mb(mt(1)))
87      value = value.and.geom_masses_set(geom2,nion2,dbl_mb(mt(1)))
88
89      if (nion3.gt.0) then
90        value = value.and.
91     >          geom_cart_set(geom3,nion3,byte_mb(tt(1)+16*nion2),
92     >                                   dbl_mb(rt(1) + 3*nion2),
93     >                                   dbl_mb(qt(1) +   nion2))
94        value = value.and.geom_masses_set(geom3,nion3,
95     >                                   dbl_mb(mt(1) +   nion2))
96      end if
97
98
99      call dcopy(nion1,0.0d0,0,dbl_mb(rt(1)),1)
100      value = value.and.geom_vel_get(geom1,dbl_mb(rt(1)))
101      value = value.and.geom_vel_set(geom2,dbl_mb(rt(1)))
102      if (nion3.gt.0)
103     > value = value.and.geom_vel_set(geom3,dbl_mb(rt(1)+3*nion2))
104
105      value = value.and.geom_rtdb_store(rtdb,geom2,'chargepspwgeometry')
106      if(nion3.gt.0) then
107         value = value.and.geom_rtdb_store(rtdb,geom3,'qmmmgeometry')
108         mmexist = .true.
109      else
110         mmexist = .false.
111      end if
112      value = value.and.geom_destroy(geom3)
113      value = value.and.geom_destroy(geom2)
114      value = value.and.geom_destroy(geom1)
115      if (.not. value)
116     >   call errquit('geometry->chargepspwgeometry write',0, GEOM_ERR)
117      value = BA_pop_stack(tt(2))
118      value = value.and.BA_pop_stack(mt(2))
119      value = value.and.BA_pop_stack(qt(2))
120      value = value.and.BA_pop_stack(rt(2))
121      if (.not. value) call errquit('popping stack',0, MA_ERR)
122
123      return
124      end
125
126*     ***************************
127*     *                         *
128*     *      parseqmmm          *
129*     *                         *
130*     ***************************
131      logical function parseqmmm(string)
132      implicit none
133      character*16 string
134
135      logical qmmm
136
137      qmmm = .false.
138      if (index(string,'^').gt.0)   qmmm = .true.
139
140      parseqmmm = qmmm
141      return
142      end
143
144
145*     ***************************
146*     *                         *
147*     *     combine_molpsp      *
148*     *                         *
149*     ***************************
150      subroutine combine_molpsp(rtdb)
151      implicit none
152#include "errquit.fh"
153      integer rtdb
154
155#include "bafdecls.fh"
156#include "btdb.fh"
157#include "geom.fh"
158
159      logical mmexist
160      common / ion_qmmm/ mmexist
161
162*     **** local variables ****
163      integer i,geom1,geom2,geom3,nion1,nion2,nion3,rt(2)
164      logical value
165      double precision rxyz(3),q
166      character*16     t
167
168
169
170*       **********************************************************
171*       **** put together ions and molecular pseudopotentials ****
172*       **********************************************************
173        value = geom_create(geom1,'geometry')
174        value = value.and.geom_create(geom2,'chargepspwgeometry')
175        if (mmexist)
176     >    value = value.and.geom_create(geom3,'qmmmgeometry')
177        value = value.and.geom_rtdb_load(rtdb,geom1,'geometry')
178        value = value.and.geom_ncent(geom1,nion1)
179        value = value.and.
180     >          geom_rtdb_load(rtdb,geom2,'chargepspwgeometry')
181        value = value.and.geom_ncent(geom2,nion2)
182        if (mmexist) then
183           value = value.and.geom_rtdb_load(rtdb,geom3,'qmmmgeometry')
184           value = value.and.geom_ncent(geom3,nion3)
185        else
186           nion3 = 0
187        end if
188        if (.not. value)
189     >   call errquit('chargepspwgeometry->geometry write 1',0,
190     &       GEOM_ERR)
191
192        value = BA_push_get(mt_dbl,(3*nion1),'rt',rt(2),rt(1))
193        if (.not. value) call errquit('out of stack memory',0, MA_ERR)
194
195        do i=1,nion2
196           value = value.and.geom_cent_get(geom2,i,t,rxyz,q)
197           value = value.and.geom_cent_set(geom1,i,t,rxyz,q)
198        end do
199        do i=1,nion3
200           value = value.and.geom_cent_get(geom3,i,t,rxyz,q)
201           value = value.and.geom_cent_set(geom1,(i+nion2),t,rxyz,q)
202        end do
203        value = value.and.geom_vel_get(geom2,dbl_mb(rt(1)))
204        if (nion3.gt.0)
205     >     value = value.and.geom_vel_get(geom3,dbl_mb(rt(1+3*nion2)))
206        value = value.and.geom_vel_set(geom1, dbl_mb(rt(1)))
207
208        value = value.and.geom_rtdb_delete(rtdb,'geometry')
209        value = value.and.geom_rtdb_delete(rtdb,'chargepspwgeometry')
210        if (nion3.gt.0)
211     >     value = value.and.geom_rtdb_delete(rtdb,'qmmmgeometry')
212        value = value.and.geom_rtdb_store(rtdb,geom1,'geometry')
213        if (mmexist)
214     >     value = value.and.geom_destroy(geom3)
215        value = value.and.geom_destroy(geom2)
216        value = value.and.geom_destroy(geom1)
217        value = value.and.BA_pop_stack(rt(2))
218        if (.not. value)
219     >   call errquit('chargepspwgeometry->geometry write 2',0,
220     &       GEOM_ERR)
221
222      return
223      end
224
225
226
227
228
229
230