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