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