1 subroutine qmmm_data_init(irtdb) 2 implicit none 3c 4#include "mafdecls.fh" 5#include "errquit.fh" 6#include "qmmm_data.fh" 7#include "qmmm_params.fh" 8#include "rtdb.fh" 9#include "global.fh" 10#include "inp.fh" 11 12 integer irtdb 13c 14 character*32 pname 15 logical qmmm_h_link 16 external qmmm_h_link 17 18 19 pname = "qmmm_data_init: " 20 21 call mm_get_tot_nqm(nqm) 22 if(ga_nodeid().eq.0) 23 > write(*,*) "number of quantum atoms",nqm 24c 25c allocate data arrays 26c -------------------- 27 call qmmm_data_alloc() 28c 29c establish persistent mapping between mm and qm atoms 30c first go pure quantum atoms and then links 31c ---------------------------------------------------- 32 call mm_get_solute_quant_ind(nqm,int_mb(i_ai)) 33 if(.not.rtdb_get(irtdb,'charge',mt_dbl,1,qcharge)) then 34 qcharge = 0.0d0 35 end if 36 37 end 38 39 subroutine qmmm_data_rdinp(irtdb) 40 implicit none 41c 42#include "mafdecls.fh" 43#include "errquit.fh" 44#include "qmmm_data.fh" 45#include "qmmm_link_data.fh" 46#include "qmmm_params.fh" 47#include "rtdb.fh" 48#include "global.fh" 49#include "inp.fh" 50#include "util.fh" 51 52 integer irtdb 53c 54 character*84 tag 55 character*32 pname 56 character*32 optimization 57 character*32 link_ecp 58 double precision eatoms 59 integer linkatm,bq_exclude 60 character*30 operation 61 character*30 region(3) 62 integer nregion 63 logical ignore 64 integer i 65 66 pname = "qmmm_data_init: " 67 68 master = ga_nodeid().eq.0 69c 70 oprint_debug = util_print("debug_info", print_debug) 71 oprint_debug = oprint_debug .and. (ga_nodeid().eq.0) 72c 73 oprint_low = util_print("low", print_low) 74c 75 eatoms = 0.0d0 76 linkatm=qmmm_hydrogen 77 bq_exclude=qmmm_no_atoms 78 link_ecp = "auto" 79c 80c initialize rtdb parameters 81 if(master) 82 > call util_print_centered(6, 83 > "QM/MM Interface Parameters",32,.true.) 84c 85c get the operation 86c ----------------- 87 tag='task:operation' 88 if (.not. rtdb_cget(irtdb, 'task:operation', 1, operation)) 89 $ operation = 'energy' 90 if (master) 91 > write(6,22) "operation",tag,operation 92c 93 tag = "qmmm:uqmatm" 94 if (.not.rtdb_get(irtdb,"qmmm:uqmatm",mt_dbl,1,eatoms)) then 95 if (.not.rtdb_put(irtdb,"qmmm:uqmatm",mt_dbl,1,eatoms)) 96 > call errquit(pname//'failed to store eatoms',0, 97 > RTDB_ERR) 98 end if 99 if(master) then 100 if(eatoms.eq.0) then 101 write(6,19) "reference energy",tag,eatoms, 102 > " <--- Warning zero value is not advisable !!!" 103 else 104 write(6,20) "reference energy",tag,eatoms 105 end if 106 end if 107 108c 109c default value for bqzone from md 110 if(.not.rtdb_get(irtdb,'md:rqmmm',mt_dbl,1,bqzone)) 111 + call errquit(pname//'failed to get bqzone from mm',0) 112c convert from nm to angstrom 113 bqzone = bqzone*cnm2ang 114c 115 tag = "qmmm:bqzone" 116 if (rtdb_get(irtdb,tag,mt_dbl,1,bqzone)) then 117 if (.not.rtdb_put(irtdb,'md:rqmmm',mt_dbl,1,bqzone/cnm2ang)) 118 > call errquit(pname//'failed to store '//tag,0, 119 > RTDB_ERR) 120 end if 121 if(master) 122 > write(6,20) "bqzone radius",tag,bqzone 123 124 125 tag ="qmmm:bq_exclude" 126 if (.not.rtdb_get(irtdb,tag,mt_int,1,bq_exclude)) then 127 if (.not.rtdb_put(irtdb,tag,mt_int,1,bq_exclude)) 128 > call errquit(pname//'failed to store bq_exclude',0, 129 > RTDB_ERR) 130 131 end if 132 if(master) then 133 if(bq_exclude.eq.qmmm_linkbond_H) then 134 write(6,22) "excluded bq's",tag, 135 > "hydrogens bonded to link atoms" 136 else if(bq_exclude.eq.qmmm_linkbond) then 137 write(6,22) "excluded bq's",tag,"all atoms bonded to links" 138 else if(bq_exclude.eq.qmmm_all_atoms) then 139 write(6,22) "excluded bq's",tag,"all" 140 else if(bq_exclude.eq.qmmm_no_atoms) then 141 write(6,22) "excluded bq's",tag,"none" 142 else 143 call errquit(pname//'invalid bq_exclude',0,RTDB_ERR) 144 end if 145 end if 146 147 tag ="qmmm:linkatm" 148 if (.not.rtdb_get(irtdb,tag,mt_int,1,linkatm)) then 149 if (.not.rtdb_put(irtdb,tag,mt_int,1,linkatm)) 150 > call errquit(pname//'failed to store'//tag,0, 151 > RTDB_ERR) 152 153 end if 154 if(master) then 155 if(linkatm.eq.qmmm_hydrogen) then 156 write(6,22) "link atom type",tag,"hydrogens" 157 else if(linkatm.eq.qmmm_halogen) then 158 write(6,22) "link atom type",tag,"halogens" 159 else 160 call errquit(pname//'invalid link atom type',0,RTDB_ERR) 161 end if 162 end if 163 164 optimization = "bfgs" 165 tag ="qmmm:optimization" 166 if (.not.rtdb_cget(irtdb,tag,1,optimization)) then 167 if (.not.rtdb_cput(irtdb,tag,1,optimization)) 168 > call errquit(pname//'failed to store'//tag,0, 169 > RTDB_ERR) 170 171 end if 172 if (master) 173 > write(6,22) "optimization method",tag,optimization 174 175c 176c region definitions 177c ------------------ 178 tag ="qmmm:region" 179 if (.not.rtdb_get(irtdb,tag(1:inp_strlen(tag))//"_n", 180 > mt_int,1,nregion)) then 181c if(operation.ne."energy" .and. 182c > operation.ne."gradient" .and. 183c > operation.ne."property" .and. 184c > operation.ne."fep" .and. 185c > operation.ne."neb" ) 186c > call errquit(pname//'cannot find '//tag,0, 187c > RTDB_ERR) 188 nregion = 0 189 end if 190 if(nregion.ne.0) then 191 if (.not.rtdb_cget(irtdb,tag,nregion,region)) 192 > call errquit(pname//"failed to get"//tag,0,RTDB_ERR) 193 end if 194 if (.not.rtdb_put(irtdb,tag(1:inp_strlen(tag))//"_n", 195 > mt_int,1,nregion)) 196 > call errquit(pname//'failed to store'//tag,0, 197 > RTDB_ERR) 198 199 if (master) then 200 tag = " " 201 do i=1,nregion 202 write(tag,'("qmmm:region ",I1)') i 203 write(6,22) tag,region(i) 204 end do 205 end if 206 207 tag ="qmmm:link_ecp" 208 if (.not.rtdb_cget(irtdb,tag,1,link_ecp)) then 209 if (.not.rtdb_cput(irtdb,tag,1,link_ecp)) 210 > call errquit(pname//'failed to store'//tag,0, 211 > RTDB_ERR) 212 213 else if (link_ecp .ne. "auto" .and. 214 > link_ecp .ne. "user" ) then 215 call errquit(pname//'unknown value'//tag,0, 216 > RTDB_ERR) 217 218 end if 219 if (master) 220 > write(6,22) "ecp on link atoms",tag,link_ecp 221 222 qmmm_interface='qm' 223 tag ="qmmm:interface" 224 ignore = rtdb_cget(irtdb,tag,1,qmmm_interface) 225 if (.not.rtdb_cget(irtdb,tag,1,qmmm_interface)) then 226 qmmm_interface = "qm" 227 end if 228 if (master) 229 > write(6,22) "interface api",tag,qmmm_interface 230 231 tag = "qmmm:bq_dynamical" 232 if (.not.rtdb_get(irtdb,tag,mt_log,1,bq_dynamical)) 233 > bq_dynamical = .false. 234 235 if (master) 236 > write(6,23) 237 238 23919 FORMAT(1X,A,T24,A,T46,F12.6,A) 24020 FORMAT(1X,A,T24,A,T46,F12.6) 24121 FORMAT(1X,A,T24,A,T46,L3) 24222 FORMAT(1X,A,T24,A,T46,A) 24323 FORMAT(1X,54("-"),//) 244 245 if (.not.rtdb_get(irtdb,'qmmm:linkatm',mt_int,1,link_atom)) 246 + call errquit('qmmm_data_init: qmmm:linkatm',link_atom, 247 & RTDB_ERR) 248 249 nqm = -1 250 251 end 252 253 subroutine qmmm_data_alloc() 254 implicit none 255c 256#include "mafdecls.fh" 257#include "errquit.fh" 258#include "qmmm_data.fh" 259c 260c 261c indexing array 262c ------------- 263 if(.not.ma_alloc_get(MT_INT, nqm, 'qmmm index array', 264 & h_ai, i_ai) ) call errquit( 265 & 'qmmm_data_alloc: unable to allocate heap space', 266 & nqm, MA_ERR) 267 call ifill(nqm,-1,int_mb(i_ai),1) 268 269 end 270 271 subroutine qmmm_data_release() 272 implicit none 273#include "mafdecls.fh" 274#include "errquit.fh" 275#include "qmmm_data.fh" 276 277 if(.not.ma_free_heap(h_ai)) 278 & call errquit('qmmm ai: Failed to deallocate heap',nqm, 279 & MA_ERR) 280 281 282 return 283 end 284 285 function qmmm_get_nqm() 286 implicit none 287#include "qmmm_data.fh" 288 289 integer qmmm_get_nqm 290 291 qmmm_get_nqm = nqm 292 293 end 294 295 function qmmm_get_i_ai() 296 implicit none 297#include "qmmm_data.fh" 298 integer qmmm_get_i_ai 299 300 qmmm_get_i_ai = i_ai 301 302 end 303 304 function qmmm_master() 305 implicit none 306#include "qmmm_data.fh" 307#include "qmmm_params.fh" 308 309 logical qmmm_master 310 311 qmmm_master = master 312 313 end 314 315 function qmmm_bq_dynamical() 316 implicit none 317#include "qmmm_data.fh" 318#include "qmmm_params.fh" 319 320 logical qmmm_bq_dynamical 321 322 qmmm_bq_dynamical = bq_dynamical 323 324 end 325 326 function qmmm_print_default(name) 327 implicit none 328#include "qmmm_data.fh" 329#include "util.fh" 330#include "global.fh" 331 332 logical qmmm_print_default 333 character*(*) name 334c 335 logical oprint 336 337 oprint = util_print(name, print_default) 338 339 qmmm_print_default = oprint .and. (ga_nodeid().eq.0) 340 341 end 342 343 function qmmm_print_debug() 344 implicit none 345#include "qmmm_data.fh" 346#include "util.fh" 347#include "global.fh" 348 349 logical qmmm_print_debug 350c 351 352c oprint = util_print("debug_info", print_debug) 353c qmmm_print_debug = oprint .and. (ga_nodeid().eq.0) 354 355 qmmm_print_debug = oprint_debug 356 end 357 358 function qmmm_print_low() 359 implicit none 360#include "qmmm_data.fh" 361#include "util.fh" 362#include "global.fh" 363 364 logical qmmm_print_low 365c 366 367 qmmm_print_low = oprint_low 368 end 369 370 function qmmm_get_bqzone() 371 implicit none 372#include "qmmm_data.fh" 373 double precision qmmm_get_bqzone 374 375 qmmm_get_bqzone = bqzone 376 377 end 378 379 function qmmm_get_interface() 380 implicit none 381#include "qmmm_data.fh" 382 character*32 qmmm_get_interface 383 384 qmmm_get_interface = qmmm_interface 385 386 end 387 388 subroutine qmmm_set_interface(a_interface) 389 implicit none 390#include "qmmm_data.fh" 391 character*(*) a_interface 392 393 qmmm_interface = a_interface 394 395 end 396 397c $Id$ 398