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