1 logical function oniom_energy(rtdb) 2* 3* $Id$ 4* 5 implicit none 6 integer rtdb 7 logical oniom 8 external oniom 9 oniom_energy = oniom(rtdb,'energy') 10 end 11 logical function oniom_gradient(rtdb) 12 implicit none 13 integer rtdb 14 logical oniom 15 external oniom 16 oniom_gradient = oniom(rtdb,'gradient') 17 end 18 logical function oniom(rtdb, task) 19 implicit none 20#include "errquit.fh" 21#include "rtdb.fh" 22#include "mafdecls.fh" 23#include "global.fh" 24#include "util.fh" 25#include "geom.fh" 26#include "inp.fh" 27#include "nwc_const.fh" 28 integer rtdb 29 character*(*) task 30c 31c Implements a two- and three-layer ONIOM model 32c 33c E(High,Real) = E(High,Model) + E(Low,Real) - E(Low,Model) 34c 35c task can currently be energy or gradient 36c 37 logical oprint, ignore, status, oprint_head 38 character*256 title 39 character*32 high_theory, low_theory, medium_theory, 40 $ gname, bname, ename 41 character*32 high_basis, low_basis, medium_basis 42 character*32 high_ecp, low_ecp, medium_ecp 43c 44 character*256 v_low_real, v_low_model, v_high_model, 45 $ v_medium_model, v_medium_inter, v_low_inter 46c 47 integer real_nat, model_nat, nat, nattest, inter_nat 48 integer geom 49 integer i, j, k, bond, nlayer, dummy(2) 50 integer maxbonds 51 parameter (maxbonds = 50) 52 integer model_bonds(2,maxbonds), model_nbonds 53 double precision model_p(maxbonds), model_charge 54 character*32 model_tags(maxbonds) 55 integer inter_bonds(2,maxbonds), inter_nbonds 56 double precision inter_p(maxbonds), inter_charge 57 character*32 inter_tags(maxbonds) 58 double precision p, real_charge 59 character*1024 high_input, medium_input, low_input 60c 61c For manipulating the geometry 62c 63 double precision coords(3,nw_max_atom), charge(nw_max_atom) 64 character*16 tags(nw_max_atom) 65c 66 double precision g(3,nw_max_atom), g_low_model(3,nw_max_atom), 67 $ g_high_model(3,nw_max_atom), g_low_real(3,nw_max_atom), 68 $ g_medium_model(3,nw_max_atom), g_medium_inter(3,nw_max_atom), 69 $ g_low_inter(3,nw_max_atom) 70c 71 double precision energy, e_low_real, e_low_model, e_high_model, 72 $ e_medium_model, e_medium_inter, e_low_inter, e2, e3 73 double precision gdummy 74c 75 character*2 symbol 76 character*8 element 77 integer atn 78c 79 integer k_actlist, l_actlist, ma_type, nactive 80c 81 logical task_energy_doit, task_gradient_doit, geom_zmtmak 82 logical usesymmol 83 integer nops 84 external task_energy_doit, task_gradient_doit, geom_zmtmak 85c 86 status = rtdb_parallel(.true.) 87 if(.not.rtdb_get(rtdb,'geom:symmol',mt_log,1, 88 $ usesymmol)) usesymmol=.false. 89 call util_print_push 90 call util_print_rtdb_load(rtdb, 'oniom') 91 call ecce_print_module_entry('oniom') 92 oprint = util_print('information', print_low) .and. 93 $ ga_nodeid().eq.0 94 oprint_head = util_print('headers', print_low) .and. 95 $ ga_nodeid().eq.0 96 if (.not. rtdb_cget(rtdb, 'title', 1, title)) 97 $ title = ' ' 98 if (.not. geom_create(geom,'geometry')) 99 $ call errquit('oniom: failed creating geometry?',0, GEOM_ERR) 100 if (.not. geom_rtdb_load(rtdb, geom, 'geometry')) 101 $ call errquit('oniom: no geometry ', 0, RTDB_ERR) 102 if (.not. geom_rtdb_store(rtdb, geom, 'oniom_geometry')) 103 $ call errquit('oniom: failed storing geometry',0, RTDB_ERR) 104 if (.not. geom_ncent(geom,real_nat)) 105 $ call errquit('oniom: number of atoms from geometry?',0, 106 & GEOM_ERR) 107 if (.not. geom_destroy(geom)) 108 $ call errquit('oniom: destroying geometry',0, GEOM_ERR) 109c 110c Save the list active list of atoms for gradient calculations 111c 112 if (task .eq. 'gradient') then 113 nactive = real_nat 114 if (rtdb_ma_get(rtdb, 'geometry:actlist', ma_type, 115 $ nactive, l_actlist)) then 116 if (.not. ma_get_index(l_actlist, k_actlist)) 117 $ call errquit('oniom: ma_get_index failed',l_actlist, 118 & MA_ERR) 119 if (.not. rtdb_put(rtdb, 'oniom:actlist', mt_int, nactive, 120 $ int_mb(k_actlist))) call errquit 121 $ ('oniom: saving actlist failed',0, RTDB_ERR) 122 if (.not. ma_free_heap(l_actlist)) call errquit 123 $ ('oniom: free of actlist?',0, MA_ERR) 124 end if 125 else 126 ignore = rtdb_delete(rtdb,'oniom:actlist') 127 end if 128c 129c Get ONIOM parameters from the database 130c 131 if (.not. rtdb_cget(rtdb, 'oniom:high:theory', 1, high_theory)) 132 $ call errquit('oniom: high_theory?', 0, RTDB_ERR) 133 if (.not. rtdb_cget(rtdb, 'oniom:medium:theory', 1,medium_theory)) 134 $ medium_theory = ' ' 135 if (.not. rtdb_cget(rtdb, 'oniom:low:theory', 1, low_theory)) 136 $ call errquit('oniom: low_theory?', 0, RTDB_ERR) 137c 138 if (.not. rtdb_cget(rtdb, 'oniom:high:basis', 1, high_basis)) 139 $ high_basis = 'ao basis' 140 if (.not. rtdb_cget(rtdb, 'oniom:medium:basis', 1, medium_basis)) 141 $ medium_basis = high_basis 142 if (.not. rtdb_cget(rtdb, 'oniom:low:basis', 1, low_basis)) 143 $ low_basis = medium_basis 144c 145 if (.not. rtdb_cget(rtdb, 'oniom:high:ecp', 1, high_ecp)) 146 $ high_ecp = ' ' 147 if (.not. rtdb_cget(rtdb, 'oniom:medium:ecp', 1, medium_ecp)) 148 $ medium_ecp = ' ' 149 if (.not. rtdb_cget(rtdb, 'oniom:low:ecp', 1, low_ecp)) 150 $ low_ecp = ' ' 151c 152 if (.not. rtdb_cget(rtdb, 'oniom:high:input', 1, high_input)) 153 $ high_input = ' ' 154 if (.not. rtdb_cget(rtdb, 'oniom:medium:input', 1, medium_input)) 155 $ medium_input = ' ' 156 if (.not. rtdb_cget(rtdb, 'oniom:low:input', 1, low_input)) 157 $ low_input = ' ' 158c 159 if (.not. rtdb_cget(rtdb, 'oniom:v_low_real', 1, 160 $ v_low_real)) 161 $ call util_file_name('lrmos', .false.,.false., v_low_real) 162 if (.not. rtdb_cget(rtdb, 'oniom:v_low_model', 1, 163 $ v_low_model)) 164 $ call util_file_name('lmmos', .false.,.false., v_low_model) 165 if (.not. rtdb_cget(rtdb, 'oniom:v_high_model', 1, 166 $ v_high_model)) 167 $ call util_file_name('hmmos', .false.,.false., v_high_model) 168 if (.not. rtdb_cget(rtdb, 'oniom:v_medium_model', 1, 169 $ v_medium_model)) 170 $ call util_file_name('mmmos', .false.,.false., v_medium_model) 171 if (.not. rtdb_cget(rtdb, 'oniom:v_medium_inter', 1, 172 $ v_medium_inter)) 173 $ call util_file_name('mimos', .false.,.false., v_medium_inter) 174 if (.not. rtdb_cget(rtdb, 'oniom:v_low_inter', 1, 175 $ v_low_inter)) 176 $ call util_file_name('limos', .false.,.false., v_low_inter) 177c 178 call util_file_name_resolve(v_low_real, .false.) 179 call util_file_name_resolve(v_low_model, .false.) 180 call util_file_name_resolve(v_high_model, .false.) 181 call util_file_name_resolve(v_medium_model, .false.) 182 call util_file_name_resolve(v_medium_inter, .false.) 183 call util_file_name_resolve(v_low_inter, .false.) 184c 185 if (.not. rtdb_get(rtdb, 'oniom:model:nat', mt_int, 1, 186 $ model_nat)) call errquit('oniom: model_nat?',0, RTDB_ERR) 187 if (.not. rtdb_get(rtdb, 'oniom:inter:nat', mt_int, 1, 188 $ inter_nat)) inter_nat = 0 189 if (.not. rtdb_get(rtdb, 'oniom:model:charge', mt_dbl, 1, 190 $ model_charge)) call errquit('oniom:model charge',0, 191 & RTDB_ERR) 192 if (.not. rtdb_get(rtdb, 'oniom:inter:charge', mt_dbl, 1, 193 $ inter_charge)) inter_charge = 0.0d0 194c 195 model_nbonds = 0 196 if (rtdb_get(rtdb,'oniom:model:nbonds',mt_int,1,model_nbonds))then 197 if (.not. rtdb_get(rtdb, 'oniom:model:bonds', mt_int, 198 $ 2*model_nbonds, model_bonds)) call errquit 199 $ ('oniom: failed to read model-real bonds?',0, 200 & RTDB_ERR) 201 if (.not. rtdb_get(rtdb, 'oniom:model:p', mt_dbl, 202 $ model_nbonds, model_p)) call errquit 203 $ ('oniom: failed to read model-real bonds factors?',0, 204 & RTDB_ERR) 205 if (.not. rtdb_cget(rtdb,'oniom:model:tags', 206 $ model_nbonds, model_tags)) call errquit 207 $ ('oniom: failed to read model-real bonds tags?',0, 208 & RTDB_ERR) 209 end if 210c 211 inter_nbonds = 0 212 if (rtdb_get(rtdb,'oniom:inter:nbonds',mt_int,1,inter_nbonds))then 213 if (.not. rtdb_get(rtdb, 'oniom:inter:bonds', mt_int, 214 $ 2*inter_nbonds, inter_bonds)) call errquit 215 $ ('oniom: failed to read inter-real bonds?',0, 216 & RTDB_ERR) 217 if (.not. rtdb_get(rtdb, 'oniom:inter:p', mt_dbl, 218 $ inter_nbonds, inter_p)) call errquit 219 $ ('oniom: failed to read inter-real bonds factors?',0, 220 & RTDB_ERR) 221 if (.not. rtdb_cget(rtdb,'oniom:inter:tags', 222 $ inter_nbonds, inter_tags)) call errquit 223 $ ('oniom: failed to read model-real bonds tags?',0, 224 & RTDB_ERR) 225 end if 226c 227c Save user indirections for the basis set and geometry 228c 229 if (.not. rtdb_cget(rtdb,'geometry',1,gname)) gname = ' ' 230 if (.not. rtdb_cget(rtdb,'ao basis',1,bname)) bname = ' ' 231 if (.not. rtdb_cget(rtdb,'ecp basis',1,ename)) ename = ' ' 232 if (.not. rtdb_get(rtdb, 'charge', mt_dbl, 1, real_charge)) 233 $ real_charge = 0.0d0 234c 235 if (model_charge .eq. -999.0) model_charge = real_charge 236 if (inter_charge .eq. -999.0) inter_charge = real_charge 237c 238 nlayer = 2 239 if (inter_nat .gt. 0) nlayer = 3 240c 241c Summarize parameters 242c 243 if (oprint) then 244 call util_print_centered(6, 'NWChem ONIOM Module', 40, .true.) 245 write(6,*) 246 write(6,*) 247 if (title .ne. ' ') then 248 call util_print_centered(6, title, 40, .false.) 249 write(6,*) 250 write(6,*) 251 end if 252c 253 write(6,10) 254 10 format(7x,'Theory',/,7x,'------') 255 if (high_ecp .eq. ' ') then 256 write(6,1) ' high', 257 $ high_theory(1:inp_strlen(high_theory)), 258 $ high_basis(1:inp_strlen(high_basis)) 259 else 260 write(6,1) ' high', 261 $ high_theory(1:inp_strlen(high_theory)), 262 $ high_basis(1:inp_strlen(high_basis)), 263 $ high_ecp(1:inp_strlen(high_ecp)) 264 end if 265 if (high_input .ne. ' ') 266 $ write(6,101) high_input(1:inp_strlen(high_input)) 267 if (medium_theory .ne. ' ') then 268 if (medium_ecp .eq. ' ') then 269 write(6,1) ' medium', 270 $ medium_theory(1:inp_strlen(medium_theory)), 271 $ medium_basis(1:inp_strlen(medium_basis)) 272 else 273 write(6,1) ' medium', 274 $ medium_theory(1:inp_strlen(medium_theory)), 275 $ medium_basis(1:inp_strlen(medium_basis)), 276 $ medium_ecp(1:inp_strlen(medium_ecp)) 277 end if 278 if (medium_input .ne. ' ') 279 $ write(6,101) medium_input(1:inp_strlen(medium_input)) 280 end if 281 if (low_ecp .eq. ' ') then 282 write(6,1) ' low', low_theory(1:inp_strlen(low_theory)), 283 $ low_basis(1:inp_strlen(low_basis)) 284 else 285 write(6,1) ' low', low_theory(1:inp_strlen(low_theory)), 286 $ low_basis(1:inp_strlen(low_basis)), 287 $ low_ecp(1:inp_strlen(low_ecp)) 288 end if 289 if (low_input .ne. ' ') 290 $ write(6,101) low_input(1:inp_strlen(low_input)) 291 1 format(6x,a7,2x,'theory=',a,2x,'basis="',a,'"',:,2x, 292 $ 'ecp="',a,'"') 293 101 format(15x,'input string="',a,'"') 294c 295 write(6,*) 296 write(6,15) 297 15 format(6x,'Vectors',/,6x,'-------') 298 if (nlayer .eq. 2) then 299 write(6,21) 'High-model', 300 $ v_high_model(1:inp_strlen(v_high_model)) 301 write(6,21) 'Low-real', 302 $ v_low_real(1:inp_strlen(v_low_real)) 303 write(6,21) 'Low-model', 304 $ v_low_model(1:inp_strlen(v_low_model)) 305 else 306 write(6,21) 'High-model', 307 $ v_high_model(1:inp_strlen(v_high_model)) 308 write(6,21) 'Medium-model', 309 $ v_medium_model(1:inp_strlen(v_medium_model)) 310 write(6,21) 'Medium-inter', 311 $ v_medium_inter(1:inp_strlen(v_medium_inter)) 312 write(6,21) 'Low-real', 313 $ v_low_real(1:inp_strlen(v_low_real)) 314 write(6,21) 'Low-inter', 315 $ v_low_inter(1:inp_strlen(v_low_inter)) 316 end if 317 21 format(8x,a12,2x,a) 318c 319 write(6,*) 320 write(6,20) 321 20 format(5x,'Geometry',/,5x,'--------') 322 write(6,2) ' real', real_nat, real_charge 323 if (inter_nat .gt. 0) write(6,2) ' inter', inter_nat, 324 $ inter_charge 325 if (inter_nbonds .gt. 0) then 326 write(6,*) ' ', 327 $ 'Bonds between intermediate and real' 328 write(6,3) (inter_bonds(1,i),inter_bonds(2,i),inter_p(i), 329 $ inter_tags(i)(1:inp_strlen(inter_tags(i))), 330 $ i=1,inter_nbonds) 331 end if 332 write(6,2) ' model', model_nat, model_charge 333 if (model_nbonds .gt. 0) then 334 write(6,*) ' Bonds between model and', 335 $ ' real/intermediate' 336 write(6,3) (model_bonds(1,i),model_bonds(2,i),model_p(i), 337 $ model_tags(i)(1:inp_strlen(model_tags(i))), 338 $ i=1,model_nbonds) 339 end if 340 2 format(6x,a7,2x,i4,' atoms with charge',f8.4) 341 3 format(20x,2i4,f6.2,2x,'"',a,'"') 342 call util_flush(6) 343 end if 344c 345c Make the model geometry resetting the symmetry & autoz info 346c 347 if (oprint) then 348 write(6,*) 349 write(6,*) ' Making model geometry ' 350 write(6,*) 351 end if 352c 353 if (.not. geom_create(geom,'geometry')) 354 $ call errquit('oniom: failed creating geometry?',0, 355 & GEOM_ERR) 356 if (.not. geom_rtdb_load(rtdb, geom, 'oniom_geometry')) 357 $ call errquit('oniom: no geometry ', 0, 358 & RTDB_ERR) 359 if (.not. geom_cart_get(geom,real_nat,tags,g,charge)) 360 $ call errquit('oniom: geom_cart_get failed', 0, GEOM_ERR) 361 call dcopy(3*real_nat, g, 1, coords, 1) 362c 363 nat = model_nat 364 do bond = 1, model_nbonds 365 i = model_bonds(1,bond) 366 j = model_bonds(2,bond) 367 p = model_p(bond) 368 nat = nat + 1 369 tags(nat) = model_tags(bond) !!!!!!!! 370 ignore = geom_tag_to_element(tags(nat), symbol, element, atn) 371 charge(nat) = atn 372 do k = 1, 3 373 coords(k,nat) = g(k,i)*(1.0d0-p) + g(k,j)*p 374 end do 375 end do 376c 377 if (.not. geom_cart_set(geom,nat,tags,coords,charge)) 378 $ call errquit('oniom: geom_cart_get failed', 0, GEOM_ERR) 379c 380 nattest = nat ! sym_nwc corrupts the input #atoms 381 call sym_nwc(geom, rtdb, nattest, .false., 1.0d0, 1d-4, 382 , nops, usesymmol) 383 if (.not. geom_ncent(geom,nattest)) call errquit 384 $ ('oniom: getting ncent',0, GEOM_ERR) 385 if (nat .ne. nattest) call errquit 386 $ ('oniom: model does not have correct symmetry', 0, GEOM_ERR) 387 ignore = geom_zmtmak(rtdb, geom, .false.) 388 if (.not. geom_rtdb_store(rtdb, geom, 'oniom model')) 389 $ call errquit('oniom: failed storing geometry',0, RTDB_ERR) 390 if (ga_nodeid().eq.0 .and. 391 $ util_print('geometry',print_default)) then 392 if (.not. geom_print(geom)) 393 $ call errquit('oniom: geom print?',0, GEOM_ERR) 394 end if 395 if (.not. geom_destroy(geom)) 396 $ call errquit('oniom: destroying geometry',0, GEOM_ERR) 397c 398c Make the intermediate geometry resetting the symmetry & autoz info 399c 400 if (inter_nat .gt. 0) then 401 if (oprint) then 402 write(6,*) 403 write(6,*) ' Making intermediate geometry ' 404 write(6,*) 405 end if 406 if (.not. geom_create(geom,'geometry')) 407 $ call errquit('oniom: failed creating geometry?',0, 408 & GEOM_ERR) 409 if (.not. geom_rtdb_load(rtdb, geom, 'oniom_geometry')) 410 $ call errquit('oniom: no geometry ', 0, RTDB_ERR) 411 if (.not. geom_cart_get(geom,real_nat,tags,g,charge)) 412 $ call errquit('oniom: geom_cart_get failed', 0, GEOM_ERR) 413 call dcopy(3*real_nat, g, 1, coords, 1) 414 nat = inter_nat 415 do bond = 1, inter_nbonds 416 i = inter_bonds(1,bond) 417 j = inter_bonds(2,bond) 418 p = inter_p(bond) 419 nat = nat + 1 420 tags(nat) = inter_tags(bond) !!!! 'H oniom' 421 charge(nat) = 1.0d0 422 do k = 1, 3 423 coords(k,nat) = g(k,i)*(1.0d0-p) + g(k,j)*p 424 end do 425 end do 426 if (.not. geom_cart_set(geom,nat,tags,coords,charge)) 427 $ call errquit('oniom: geom_cart_get failed', 0, GEOM_ERR) 428 nattest = nat ! sym_nwc corrupts the input #atoms 429 call sym_nwc(geom, rtdb, nattest, .false., 1.0d0, 1d-4, 430 , nops, usesymmol) 431 if (.not. geom_ncent(geom,nattest)) call errquit 432 $ ('oniom: getting ncent',0, GEOM_ERR) 433 if (nat .ne. nattest) call errquit 434 $ ('oniom: inter does not have correct symmetry', 0, 435 & GEOM_ERR) 436 ignore = geom_zmtmak(rtdb, geom, .false.) 437 if (.not. geom_rtdb_store(rtdb, geom, 'oniom inter')) 438 $ call errquit('oniom: failed storing geometry',0, RTDB_ERR) 439 if (ga_nodeid().eq.0 .and. 440 $ util_print('geometry',print_default)) then 441 if (.not. geom_print(geom)) 442 $ call errquit('oniom: geom print?',0, GEOM_ERR) 443 end if 444 if (.not. geom_destroy(geom)) 445 $ call errquit('oniom: destroying geometry',0, 446 & GEOM_ERR) 447 end if 448c 449c Two-layer model 450c E(High,Real) = E(High,Model) + E(Low,Real) - E(Low,Model) 451c 452c Three-layer model 453c E(High,Real) = E(High,Model) - E(Medium,Model) 454c . + E(Medium,Inter) + E(Real,Low) - E(Low,Inter) 455c 456c To avoid recursive calls thru Fortran, call the task*doit layer. 457c 458 oniom = .false. 459 call dfill(3*nw_max_atom, 0.0d0, g, 1) 460 call dfill(3*nw_max_atom, 0.0d0, g_low_model, 1) 461 call dfill(3*nw_max_atom, 0.0d0, g_high_model, 1) 462 call dfill(3*nw_max_atom, 0.0d0, g_low_real, 1) 463 call dfill(3*nw_max_atom, 0.0d0, g_medium_model, 1) 464 call dfill(3*nw_max_atom, 0.0d0, g_medium_inter, 1) 465 energy = 0.0d0 466 e_low_model = 0.0d0 467 e_high_model = 0.0d0 468 e_low_real = 0.0d0 469 e_medium_model = 0.0d0 470 e_medium_inter = 0.0d0 471c 472c Low + Real 473c 474 if (oprint_head) 475 $ call util_print_centered(6, 'ONIOM LOW+REAL', 20, .true.) 476 call oniom_param(rtdb, ' ', low_basis, low_ecp, real_charge, 477 $ low_input) 478 call oniom_manage_mos(rtdb,low_theory,v_low_real) 479 if (task .eq. 'energy') then 480 status = task_energy_doit(rtdb, low_theory, e_low_real) 481 else 482 call oniom_munge_actlist(rtdb,real_nat,0,dummy) 483 status = task_gradient_doit(rtdb, low_theory, e_low_real, 484 $ g_low_real) 485 end if 486 if (.not. status) 487 $ call errquit('oniom: low theory + real geometry failed',0, 488 & GEOM_ERR) 489c 490c Low + Model 491c 492 if (oprint_head) 493 $ call util_print_centered(6, 'ONIOM LOW+MODEL', 20, .true.) 494 call oniom_param(rtdb, 'oniom model',low_basis, low_ecp, 495 $ model_charge, low_input) 496 call oniom_manage_mos(rtdb,low_theory,v_low_model) 497 if (task .eq. 'energy') then 498 status = task_energy_doit(rtdb, low_theory, e_low_model) 499 else 500 call oniom_munge_actlist(rtdb,model_nat, 501 $ model_nbonds,model_bonds) 502 status = task_gradient_doit(rtdb, low_theory, e_low_model, 503 $ g_low_model) 504 end if 505 if (.not. status) 506 $ call errquit('oniom:low theory + model geometry failed',0, 507 & GEOM_ERR) 508c 509c High + Model 510c 511 if (oprint_head) 512 $ call util_print_centered(6, 'ONIOM HIGH+MODEL', 20, .true.) 513 call oniom_param(rtdb, 'oniom model', 514 $ high_basis, high_ecp, model_charge, high_input) 515 call oniom_manage_mos(rtdb,high_theory, v_high_model) 516 if (task .eq. 'energy') then 517 status = task_energy_doit(rtdb, high_theory, e_high_model) 518 else 519 call oniom_munge_actlist(rtdb,model_nat, 520 $ model_nbonds,model_bonds) 521 status = task_gradient_doit(rtdb, high_theory, e_high_model, 522 $ g_high_model) 523 end if 524 if (.not. status) 525 $ call errquit('oniom: high theory + model geometry failed',0, 526 & GEOM_ERR) 527c 528 if (nlayer .eq. 3) then 529c 530c Low + Inter 531c 532 if (oprint_head) 533 $ call util_print_centered(6, 'ONIOM LOW+INTER', 20, .true.) 534 call oniom_param(rtdb, 'oniom inter', 535 $ low_basis, low_ecp, inter_charge, low_input) 536 call oniom_manage_mos(rtdb, low_theory, v_low_inter) 537 if (task .eq. 'energy') then 538 status = task_energy_doit(rtdb, low_theory, e_low_inter) 539 else 540 call oniom_munge_actlist(rtdb,inter_nat, 541 $ inter_nbonds,inter_bonds) 542 status = task_gradient_doit(rtdb, low_theory, e_low_inter, 543 $ g_low_inter) 544 end if 545 if (.not. status) 546 $ call errquit('oniom:low theory + model geometry failed',0, 547 & GEOM_ERR) 548c 549c Medium + Inter 550c 551 if (oprint_head) call util_print_centered 552 $ (6, 'ONIOM MEDIUM+INTER', 20, .true.) 553 call oniom_param(rtdb, 'oniom inter', 554 $ medium_basis, medium_ecp, inter_charge, medium_input) 555 call oniom_manage_mos(rtdb,medium_theory,v_medium_inter) 556 if (task .eq. 'energy') then 557 status = task_energy_doit(rtdb,medium_theory,e_medium_inter) 558 else 559 call oniom_munge_actlist(rtdb,inter_nat, 560 $ inter_nbonds,inter_bonds) 561 status = task_gradient_doit(rtdb, medium_theory, 562 $ e_medium_inter, g_medium_inter) 563 end if 564 if (.not. status) 565 $ call errquit('oniom: medium theory + '// 566 $ 'inter geometry failed',0, GEOM_ERR) 567c 568c Medium + Model 569c 570 if (oprint_head) call util_print_centered 571 $ (6, 'ONIOM MEDIUM+MODEL', 20, .true.) 572 call oniom_param(rtdb, 'oniom model', 573 $ medium_basis, medium_ecp, model_charge, medium_input) 574 call oniom_manage_mos(rtdb,medium_theory,v_medium_model) 575 if (task .eq. 'energy') then 576 status = task_energy_doit(rtdb,medium_theory,e_medium_model) 577 else 578 call oniom_munge_actlist(rtdb,model_nat, 579 $ model_nbonds,model_bonds) 580 status = task_gradient_doit(rtdb, medium_theory, 581 $ e_medium_model, g_medium_model) 582 end if 583 if (.not. status) 584 $ call errquit('oniom: medium theory + '// 585 $ 'model geometry failed',0, GEOM_ERR) 586 end if 587c 588c Assemble the ONIOM energy and gradient 589c 590 e2 = e_high_model + e_low_real - e_low_model 591 energy = e2 592 if (nlayer .eq. 3) then 593 e3 = e_high_model - e_medium_model 594 $ + e_medium_inter + e_low_real - e_low_inter 595 energy = e3 596 end if 597c 598 if (.not. rtdb_put(rtdb,'oniom:energy',mt_dbl,1,energy)) 599 $ call errquit('oniom: failed storing energy', 0, RTDB_ERR) 600c 601 if (task .eq. 'gradient') then 602 if (nlayer .eq. 2) then 603 call dcopy(3*real_nat, g_low_real, 1, g, 1) 604 do i = 1, model_nat 605 do k = 1, 3 606 g(k,i) = g(k,i) + g_high_model(k,i) - g_low_model(k,i) 607 end do 608 end do 609 nat = model_nat 610 do bond = 1, model_nbonds 611 i = model_bonds(1,bond) 612 j = model_bonds(2,bond) 613 p = model_p(bond) 614 nat = nat + 1 615 do k = 1, 3 616 gdummy = g_high_model(k,nat) - g_low_model(k,nat) 617 g(k,i) = g(k,i) + gdummy*(1.0d0-p) 618 g(k,j) = g(k,j) + gdummy*p 619 end do 620 end do 621 else 622 call dcopy(3*real_nat, g_low_real, 1, g, 1) 623 do i = 1, model_nat 624 do k = 1, 3 625 g(k,i) = g(k,i) + 626 $ g_high_model(k,i) - g_medium_model(k,i) 627 end do 628 end do 629 nat = model_nat 630 do bond = 1, model_nbonds 631 i = model_bonds(1,bond) 632 j = model_bonds(2,bond) 633 p = model_p(bond) 634 nat = nat + 1 635 do k = 1, 3 636 gdummy = g_high_model(k,nat) - g_medium_model(k,nat) 637 g(k,i) = g(k,i) + gdummy*(1.0d0-p) 638 g(k,j) = g(k,j) + gdummy*p 639 end do 640 end do 641 do i = 1, inter_nat 642 do k = 1, 3 643 g(k,i) = g(k,i) + 644 $ g_medium_inter(k,i) - g_low_inter(k,i) 645 end do 646 end do 647 nat = inter_nat 648 do bond = 1, inter_nbonds 649 i = inter_bonds(1,bond) 650 j = inter_bonds(2,bond) 651 p = inter_p(bond) 652 nat = nat + 1 653 do k = 1, 3 654 gdummy = g_medium_inter(k,nat) - g_low_inter(k,nat) 655 g(k,i) = g(k,i) + gdummy*(1.0d0-p) 656 g(k,j) = g(k,j) + gdummy*p 657 end do 658 end do 659 end if 660c 661c Enforce active atom constraints on the gradient and also restore 662c the full list of active atoms from the copy made earlier. 663c 664 if (rtdb_ma_get(rtdb, 'oniom:actlist', ma_type, 665 $ nactive, l_actlist)) then 666 if (.not. ma_get_index(l_actlist, k_actlist)) 667 $ call errquit('oniom: ma_get_index failed',l_actlist, 668 & MA_ERR) 669 do i = 1, real_nat 670 do j = 1, nactive 671 if (int_mb(k_actlist+j-1) .eq. i) goto 444 672 end do 673 do k = 1, 3 674 g(k,i) = 0.0d0 675 end do 676 444 continue 677 end do 678 if (.not. rtdb_put(rtdb, 'geometry:actlist', mt_int, 679 $ nactive, int_mb(k_actlist))) call errquit 680 $ ('oniom: saving actlist failed',0, 681 & RTDB_ERR) 682 if (.not. ma_free_heap(l_actlist)) call errquit 683 $ ('oniom: free of actlist?',0, RTDB_ERR) 684 end if 685c 686 if (.not. rtdb_put(rtdb,'oniom:gradient',mt_dbl,3*real_nat,g)) 687 $ call errquit('oniom: failed storing gradient', 0, 688 & RTDB_ERR) 689 end if 690c 691c Reset the geometry and restore any user indirection of basis 692c or geometry. 693c 694 ignore = rtdb_delete(rtdb, 'geometry') 695 ignore = rtdb_delete(rtdb, 'ao basis') 696 ignore = rtdb_delete(rtdb, 'ecp basis') 697 ignore = rtdb_delete(rtdb, 'charge') 698 if (gname .ne. ' ') then 699 if (.not. rtdb_cput(rtdb,'geometry',1,gname)) 700 $ call errquit('oniom:failed resetting user geom name',0, 701 & RTDB_ERR) 702 end if 703 if (bname .ne. ' ') then 704 if (.not. rtdb_cput(rtdb,'ao basis',1,bname)) call errquit 705 $ ('oniom:failed resetting user basis name',0, RTDB_ERR) 706 end if 707 if (ename .ne. ' ') then 708 if (.not. rtdb_cput(rtdb,'ecp basis',1,ename)) call errquit 709 $ ('oniom:failed resetting user ecp name',0, 710 & RTDB_ERR) 711 end if 712 if (real_charge .ne. -999.0d0) then 713 if (.not. rtdb_put(rtdb, 'charge', mt_dbl, 1, real_charge)) 714 $ call errquit('oniom: resetting real charge?', 0, 715 & RTDB_ERR) 716 end if 717c 718c Don't leave MO vectors info lying around in the database. 719c 720 call oniom_unmanage_mos(rtdb) 721c 722c Get the geometry in preparation for printing 723c 724 if (.not. geom_create(geom,'geometry')) 725 $ call errquit('oniom: failed creating geometry?',0, GEOM_ERR) 726 if (.not. geom_rtdb_load(rtdb, geom, 'oniom_geometry')) 727 $ call errquit('oniom: no geometry ', 0, RTDB_ERR) 728 if (.not. geom_cart_get(geom,real_nat,tags,coords,charge)) 729 $ call errquit('oniom: geom_cart_get failed', 0, GEOM_ERR) 730 if (.not. geom_destroy(geom)) 731 $ call errquit('oniom: geom destroy?',0, GEOM_ERR) 732c 733 if (oprint) then 734 call util_print_centered(6, 'ONIOM summary', 20, .true.) 735 write(6,*) 736 write(6,5) e_low_real, e_high_model, e_low_model, e2 737 5 format( 738 $ 8x, '+ LOW + REAL ', f20.12,/, 739 $ 8x, '+ HIGH + MODEL ', f20.12,/, 740 $ 8x, '- LOW + MODEL ', f20.12,/, 741 $ 8x, ' ------------ ',/, 742 $ 8x, '= ONIOM2 ', f20.12,/) 743 if (nlayer .eq. 3) then 744 write(6,55) e_low_real, e_high_model, e_medium_model, 745 $ e_medium_inter, e_low_inter, e3 746 55 format( 747 $ 8x, '+ LOW + REAL ', f20.12,/, 748 $ 8x, '+ HIGH + MODEL ', f20.12,/, 749 $ 8x, '- MED + MODEL ', f20.12,/, 750 $ 8x, '+ MED + INTER ', f20.12,/, 751 $ 8x, '- LOW + INTER ', f20.12,/, 752 $ 8x, ' ------------ ',/, 753 $ 8x, '= ONIOM3 ', f20.12,/) 754 end if 755 if (task .eq. 'gradient') then 756 write(6,1000) 'ONIOM', 757 $ 'x','y','z','x','y','z' 758 do 30, i=1, real_nat 759 write(6,2000) i, tags(i),(coords(j,i),j=1,3), 760 $ (g(j,i),j=1,3) 761 30 continue 762 write(6,*) 763 1000 format(/,/,25X,A,' ENERGY GRADIENTS',/,/,4X,'atom',15X, 764 $ 'coordinates', 765 $ 24X,'gradient',/,6X,2(1X,(3(10X,A1)))) 766 2000 format(1X,I3,1X,A4,2(1X,3(1X,F10.6))) 767 end if 768 call util_flush(6) 769 end if 770c 771 if (.not. rtdb_cput(rtdb, 'task:theory', 1, 'oniom')) 772 $ call errquit('oniom: re-setting theory?',0, 773 & RTDB_ERR) 774c 775 oniom = .true. 776 call util_print_pop 777 call ecce_print_module_exit('oniom', 'ok') 778c 779 end 780 subroutine oniom_input(rtdb) 781 implicit none 782#include "errquit.fh" 783#include "inp.fh" 784#include "mafdecls.fh" 785#include "rtdb.fh" 786#include "geom.fh" 787c 788 integer rtdb 789c 790 logical status 791 character*32 test 792 integer maxbonds 793 parameter (maxbonds = 50) 794c 795 integer model_bonds(2,maxbonds), model_nbonds, model_nat 796 double precision model_p(maxbonds), model_charge 797 character*32 model_tags(maxbonds) 798c 799 integer inter_bonds(2,maxbonds), inter_nbonds, inter_nat 800 double precision inter_p(50), inter_charge 801 character*32 inter_tags(maxbonds) 802c 803 character*32 high_theory, high_basis, high_ecp 804 character*32 medium_theory, medium_basis, medium_ecp 805 character*32 low_theory, low_basis, low_ecp 806 character*1024 high_input, medium_input, low_input 807c 808 character*256 v_low_real, v_low_model, v_high_model, 809 $ v_medium_model, v_medium_inter, v_low_inter 810c 811 integer bond, junk, i, j 812 character*2 symbol 813 character*8 element 814 integer atn 815c 816 high_theory = ' ' 817 high_basis = ' ' 818 high_ecp = ' ' 819 high_input = ' ' 820 medium_theory = ' ' 821 medium_basis = ' ' 822 medium_ecp = ' ' 823 medium_input = ' ' 824 low_theory = ' ' 825 low_basis = ' ' 826 low_ecp = ' ' 827 low_input = ' ' 828 model_nat = 0 829 model_nbonds = 0 830 model_charge = -999.0 ! Used to indicate no charge 831 inter_nat = 0 832 inter_nbonds = 0 833 inter_charge = -999.0 834c 835 v_low_real = ' ' 836 v_low_model = ' ' 837 v_high_model = ' ' 838 v_medium_model = ' ' 839 v_medium_inter = ' ' 840 v_low_inter = ' ' 841c 842c While input is available 843c 844 10 if (.not. inp_read()) 845 $ call errquit('oniom_input: premature EOF',0, INPUT_ERR) 846 if (.not. inp_a(test)) 847 $ call errquit('oniom_input: failed to read keyword',0, 848 & INPUT_ERR) 849c 850 if (inp_compare(.false.,'high',test)) then 851c 852 call oniom_theory_input(high_theory,high_basis,high_ecp, 853 $ high_input) 854c 855 else if (inp_compare(.false.,'medium',test)) then 856c 857 call oniom_theory_input(medium_theory,medium_basis,medium_ecp, 858 $ medium_input) 859c 860 else if (inp_compare(.false.,'low',test)) then 861c 862 call oniom_theory_input(low_theory,low_basis,low_ecp, 863 $ low_input) 864c 865 else if (inp_compare(.false.,'print', test) .or. 866 $ inp_compare(.false.,'noprint', test)) then 867 call util_print_input(rtdb, 'oniom') 868c 869 else if (inp_compare(.false.,'model',test)) then 870c 871c model natom [charge q] [i1 j1 p1 [tag1] ...] 872c 873 if (.not. inp_i(model_nat)) 874 $ call errquit('oniom_input:model natoms [i1 j1 p1 ...]',0, 875 & INPUT_ERR) 876 if (inp_a(test)) then 877 if (inp_compare(.false.,'charge',test)) then 878 if (.not. inp_f(model_charge)) call errquit 879 $ ('oniom_input: reading charge', 0, INPUT_ERR) 880 else 881 call inp_prev_field() 882 end if 883 end if 884 model_nbonds = 0 885 do bond = 1, maxbonds 886 if (.not. inp_i(model_bonds(1,bond))) goto 501 887 if (.not. inp_i(model_bonds(2,bond))) call errquit 888 $ ('oniom_input: reading second atom for bond',bond, 889 & INPUT_ERR) 890 if (.not. inp_f(model_p(bond))) call errquit 891 $ ('oniom_input: reading scale factor for bond',bond, 892 & INPUT_ERR) 893 model_tags(bond) = 'H oniom' 894 if (inp_cur_field() .lt. inp_n_field()) then 895 if (inp_i(junk)) then 896 call inp_prev_field() 897 else if (.not. inp_a(model_tags(bond))) then 898 call errquit('oniom_input: reading tag for bond',bond, 899 & INPUT_ERR) 900 end if 901 end if 902 if (.not. geom_tag_to_element(model_tags(bond), symbol, 903 $ element, atn)) then 904 if (symbol.ne.'bq') call errquit 905 $ ('oniom_input: tag must be bq* or element*',0, 906 & INPUT_ERR) 907 end if 908 model_nbonds = model_nbonds + 1 909 end do 910 501 continue 911 else if (inp_compare(.false.,'inter',test)) then 912c 913c inter natom [i1 j1 p1 [tag1] ...] 914c 915 if (.not. inp_i(inter_nat)) 916 $ call errquit('oniom_input:inter natoms [i1 j1 p1 ...]',0, 917 & INPUT_ERR) 918 if (inp_a(test)) then 919 if (inp_compare(.false.,'charge',test)) then 920 if (.not. inp_f(inter_charge)) call errquit 921 $ ('oniom_input: reading charge', 0, INPUT_ERR) 922 else 923 call inp_prev_field() 924 end if 925 end if 926 inter_nbonds = 0 927 do bond = 1, maxbonds 928 if (.not. inp_i(inter_bonds(1,bond))) goto 503 929 if (.not. inp_i(inter_bonds(2,bond))) call errquit 930 $ ('oniom_input: reading second atom for bond',bond, 931 & INPUT_ERR) 932 if (.not. inp_f(inter_p(bond))) call errquit 933 $ ('oniom_input: reading scale factor for bond',bond, 934 & INPUT_ERR) 935 inter_tags(bond) = 'H oniom' 936 if (inp_cur_field() .lt. inp_n_field()) then 937 if (inp_i(junk)) then 938 call inp_prev_field() 939 else if (.not. inp_a(inter_tags(bond))) then 940 call errquit('oniom_input: reading tag for bond',bond, 941 & INPUT_ERR) 942 end if 943 end if 944 if (.not. geom_tag_to_element(inter_tags(bond), symbol, 945 $ element, atn)) then 946 if (symbol.ne.'bq') call errquit 947 $ ('oniom_input: tag must be bq* or element*',0, 948 & INPUT_ERR) 949 end if 950 inter_nbonds = inter_nbonds + 1 951 end do 952 503 continue 953c 954 else if (inp_compare(.false.,'vectors',test)) then 955 111 if (inp_a(test)) then 956 if (inp_compare(.false.,'low-real',test)) then 957 if (.not. inp_a(v_low_real)) call errquit 958 $ ('oniom_input: error reading low_real vectors', 0, 959 & INPUT_ERR) 960 else if (inp_compare(.false.,'low-model',test)) then 961 if (.not. inp_a(v_low_model)) call errquit 962 $ ('oniom_input: error reading low_model vectors', 0, 963 & INPUT_ERR) 964 else if (inp_compare(.false.,'high-model',test)) then 965 if (.not. inp_a(v_high_model)) call errquit 966 $ ('oniom_input: error reading high_model vectors', 0, 967 & INPUT_ERR) 968 else if (inp_compare(.false.,'medium-model',test)) then 969 if (.not. inp_a(v_medium_model)) call errquit 970 $ ('oniom_input:error reading medium_model vectors',0, 971 & INPUT_ERR) 972 else if (inp_compare(.false.,'medium-inter',test)) then 973 if (.not. inp_a(v_medium_inter)) call errquit 974 $ ('oniom_input:error reading medium_inter vectors',0, 975 & INPUT_ERR) 976 else if (inp_compare(.false.,'low-inter',test)) then 977 if (.not. inp_a(v_low_inter)) call errquit 978 $ ('oniom_input: error reading low_inter vectors', 0, 979 & INPUT_ERR) 980 else 981 call errquit('oniom_input:unknown keyword for vectors',0, 982 & INPUT_ERR) 983 end if 984 goto 111 985 end if 986 else if (.not. inp_compare(.false.,'end',test)) then 987 call errquit('oniom_input: unknown directive', 0, INPUT_ERR) 988 end if 989 if (.not. inp_compare(.false.,'end',test)) goto 10 990c 991c Sanity tests 992c 993 if (high_theory.eq.' ' .or. low_theory.eq.' ') call errquit 994 $ ('oniom_input: must specify both high and low theory',0, 995 & INPUT_ERR) 996 if (model_nat .le. 0) call errquit 997 $ ('oniom_input: must specify no. of model atoms',0, INPUT_ERR) 998 if ( (medium_theory.ne.' ' .and. inter_nat.eq.0) .or. 999 $ (medium_theory.eq.' ' .and. inter_nat.gt.0) ) call errquit 1000 $ ('oniom_input: 3-layer model requires both medium theory'// 1001 $ ' and intermediate layer',0, INPUT_ERR) 1002c 1003c Put the bond info into the right order 1=model, 2=real and 1004c also do a little checking. 1005c 1006 do bond = 1, model_nbonds 1007 i = min(model_bonds(1,bond),model_bonds(2,bond)) 1008 j = max(model_bonds(1,bond),model_bonds(2,bond)) 1009 if (i.le.0 .or. i.gt.model_nat .or. j.le.model_nat) 1010 $ call errquit 1011 $ ('oniom_input: bad atom info for model bond number', bond, 1012 & INPUT_ERR) 1013 model_bonds(1,bond) = i 1014 model_bonds(2,bond) = j 1015 end do 1016 do bond = 1, inter_nbonds 1017 i = min(inter_bonds(1,bond),inter_bonds(2,bond)) 1018 j = max(inter_bonds(1,bond),inter_bonds(2,bond)) 1019 if (i.le.0 .or. i.gt.inter_nat .or. j.le.inter_nat) 1020 $ call errquit 1021 $ ('oniom_input: bad atom info for inter bond number', bond, 1022 & INPUT_ERR) 1023 inter_bonds(1,bond) = i 1024 inter_bonds(2,bond) = j 1025 end do 1026c 1027c Store info to the database 1028c 1029 status = rtdb_cput(rtdb, 'oniom:high:theory', 1, high_theory) 1030 status = rtdb_cput(rtdb, 'oniom:low:theory', 1, low_theory) 1031 $ .and. status 1032 if (medium_theory .ne. ' ') 1033 $ status = rtdb_cput(rtdb, 'oniom:medium:theory', 1, 1034 $ medium_theory) .and. status 1035c 1036 if (high_basis .ne. ' ') 1037 $ status = rtdb_cput(rtdb, 'oniom:high:basis', 1, high_basis) 1038 $ .and. status 1039 if (medium_basis .ne. ' ') 1040 $ status = rtdb_cput(rtdb, 'oniom:medium:basis', 1, 1041 $ medium_basis).and. status 1042 if (low_basis .ne. ' ') 1043 $ status = rtdb_cput(rtdb, 'oniom:low:basis', 1, low_basis) 1044 $ .and. status 1045c 1046 if (high_ecp .ne. ' ') 1047 $ status = rtdb_cput(rtdb, 'oniom:high:ecp', 1, high_ecp) 1048 $ .and. status 1049 if (medium_ecp .ne. ' ') 1050 $ status = rtdb_cput(rtdb, 'oniom:medium:ecp', 1, 1051 $ medium_ecp).and. status 1052 if (low_ecp .ne. ' ') 1053 $ status = rtdb_cput(rtdb, 'oniom:low:ecp', 1, low_ecp) 1054 $ .and. status 1055c 1056 if (high_input .ne. ' ') 1057 $ status = rtdb_cput(rtdb, 'oniom:high:input', 1, high_input) 1058 $ .and. status 1059 if (medium_input .ne. ' ') 1060 $ status = rtdb_cput(rtdb, 'oniom:medium:input', 1, 1061 $ medium_input).and. status 1062 if (low_input .ne. ' ') 1063 $ status = rtdb_cput(rtdb, 'oniom:low:input', 1, low_input) 1064 $ .and. status 1065c 1066 if (v_low_real .ne. ' ') status = status .and. 1067 $ rtdb_cput(rtdb, 'oniom:v_low_real', 1, v_low_real) 1068 if (v_low_model .ne. ' ') status = status .and. 1069 $ rtdb_cput(rtdb, 'oniom:v_low_model', 1, v_low_model) 1070 if (v_high_model .ne. ' ') status = status .and. 1071 $ rtdb_cput(rtdb, 'oniom:v_high_model', 1, v_high_model) 1072 if (v_medium_model .ne. ' ') status = status .and. 1073 $ rtdb_cput(rtdb, 'oniom:v_medium_model', 1, v_medium_model) 1074 if (v_medium_inter .ne. ' ') status = status .and. 1075 $ rtdb_cput(rtdb, 'oniom:v_medium_inter', 1, v_medium_inter) 1076 if (v_low_inter .ne. ' ') status = status .and. 1077 $ rtdb_cput(rtdb, 'oniom:v_low_inter', 1, v_low_inter) 1078c 1079 status = rtdb_put(rtdb, 'oniom:model:nat', mt_int, 1, model_nat) 1080 $ .and. status 1081 status = rtdb_put(rtdb, 'oniom:model:charge', mt_dbl, 1, 1082 $ model_charge) .and. status 1083 if (model_nbonds .gt. 0) then 1084 status = rtdb_put(rtdb, 'oniom:model:nbonds', mt_int, 1, 1085 $ model_nbonds) .and. status 1086 status = rtdb_put(rtdb, 'oniom:model:bonds', mt_int, 1087 $ 2*model_nbonds, model_bonds) .and. status 1088 status = rtdb_put(rtdb, 'oniom:model:p', mt_dbl, 1089 $ model_nbonds, model_p) .and. status 1090 status = rtdb_cput(rtdb,'oniom:model:tags', 1091 $ model_nbonds, model_tags) .and. status 1092 end if 1093c 1094 if (inter_nat .gt. 0) then 1095 status = rtdb_put(rtdb, 'oniom:inter:nat', mt_int, 1, 1096 $ inter_nat) .and. status 1097 status = rtdb_put(rtdb, 'oniom:inter:charge', mt_dbl, 1, 1098 $ inter_charge) .and. status 1099 if (inter_nbonds .gt. 0) then 1100 status = rtdb_put(rtdb, 'oniom:inter:nbonds', mt_int, 1, 1101 $ inter_nbonds) .and. status 1102 status = rtdb_put(rtdb, 'oniom:inter:bonds', mt_int, 1103 $ 2*inter_nbonds, inter_bonds) .and. status 1104 status = rtdb_put(rtdb, 'oniom:inter:p', mt_dbl, 1105 $ inter_nbonds, inter_p) .and. status 1106 status = rtdb_cput(rtdb,'oniom:inter:tags', 1107 $ inter_nbonds, inter_tags) .and. status 1108 end if 1109 end if 1110c 1111 if (.not. status) call errquit 1112 $ ('oniom_input: error putting data into RTDB',0, RTDB_ERR) 1113c 1114 end 1115 subroutine oniom_unmanage_mos(rtdb) 1116#include "rtdb.fh" 1117 integer rtdb 1118c 1119c Delete the vectors input/output keywords from the 1120c popular (scf,dft,mcscf) places so that there is 1121c no confusion or overwriting down the road 1122c 1123 logical ignore 1124c 1125 ignore = rtdb_delete(rtdb, 'scf:input vectors') 1126 ignore = rtdb_delete(rtdb, 'dft:input vectors') 1127 ignore = rtdb_delete(rtdb, 'pspw:input vectors') 1128 ignore = rtdb_delete(rtdb, 'band:input vectors') 1129 ignore = rtdb_delete(rtdb, 'mcscf:input vectors') 1130 ignore = rtdb_delete(rtdb, 'scf:output vectors') 1131 ignore = rtdb_delete(rtdb, 'dft:output vectors') 1132 ignore = rtdb_delete(rtdb, 'pspw:output vectors') 1133 ignore = rtdb_delete(rtdb, 'band:output vectors') 1134 ignore = rtdb_delete(rtdb, 'mcscf:output vectors') 1135c 1136 end 1137 subroutine oniom_manage_mos(rtdb,theory,movecs) 1138 implicit none 1139#include "errquit.fh" 1140#include "rtdb.fh" 1141#include "global.fh" 1142#include "inp.fh" 1143 integer rtdb 1144 character*(*) theory 1145 character*(*) movecs 1146c 1147c Defeat the convergence tests on the MO methods which 1148c are confused by some ONIOM calculations and tell the 1149c MO codes to use the correct MO vectors. 1150c 1151 logical ignore 1152 logical oexists 1153 logical do_mcscf 1154 character*32 keyin, keyout 1155c 1156 do_mcscf = .false. 1157 if (ga_nodeid() .eq. 0) then 1158 keyin = 'scf:input vectors' 1159 keyout = 'scf:output vectors' 1160 if (inp_compare(.false.,'dft',theory)) then 1161 keyin = 'dft:input vectors' 1162 keyout = 'dft:output vectors' 1163 else if (inp_compare(.false.,'pspw',theory)) then 1164 keyin = 'pspw:input vectors' 1165 keyout = 'pspw:output vectors' 1166 else if (inp_compare(.false.,'band',theory)) then 1167 keyin = 'band:input vectors' 1168 keyout = 'band:output vectors' 1169 else if (inp_compare(.false.,'mcscf',theory)) then 1170 keyin = 'mcscf:input vectors' 1171 keyout = 'mcscf:output vectors' 1172 do_mcscf = .true. 1173 else if (inp_compare(.false.,'selci',theory)) then 1174 keyin = 'mcscf:input vectors' 1175 keyout = 'mcscf:output vectors' 1176 do_mcscf = .true. 1177 end if 1178 inquire(file=movecs,exist=oexists) 1179c 1180 ignore = rtdb_parallel(.false.) 1181 ignore = rtdb_delete(rtdb, 'dft:converged') ! DO WE REALLY NEED TO DO THIS? 1182 ignore = rtdb_delete(rtdb, 'scf:converged') 1183 ignore = rtdb_delete(rtdb, 'mcscf:converged') 1184 ignore = rtdb_delete(rtdb, keyin) 1185 ignore = rtdb_delete(rtdb, keyout) 1186 if (oexists) then 1187 if (.not. rtdb_cput(rtdb,keyin,1,movecs)) call errquit 1188 $ ('oniom_manage_mos: failed settting input vectors',0, 1189 & RTDB_ERR) 1190 else if (do_mcscf) then 1191c we need to prevent the MCSCF code picking up random SCF vectors 1192 if (.not. rtdb_cput(rtdb,keyin,1,'atomic')) call errquit 1193 $ ('oniom_manage_mos: failed settting input vectors',0, 1194 & RTDB_ERR) 1195 end if 1196 if (.not. rtdb_cput(rtdb,keyout,1,movecs)) call errquit 1197 $ ('oniom_manage_mos: failed settting output vectors',0, 1198 & RTDB_ERR) 1199 ignore = rtdb_parallel(.true.) 1200 end if 1201c 1202 call ga_sync() 1203c 1204 end 1205 subroutine oniom_munge_actlist(rtdb,nat,nbonds,bonds) 1206 implicit none 1207#include "errquit.fh" 1208#include "rtdb.fh" 1209#include "mafdecls.fh" 1210#include "nwc_const.fh" 1211 integer rtdb 1212 integer nat, nbonds, bonds(2,*) 1213c 1214c For constrained geometry optimizations and numerical frequencies 1215c the list of active atoms for computing the gradient may be set. 1216c For ONIOM subcalculations we need to prune the list of atoms 1217c to those actually active in the calculation (otherwise some of 1218c the other modules rightfully complain) and if either of the 1219c atoms in a bond are active we need to ensure that the link 1220c atom is also marked active. 1221c 1222 integer l_actlist, k_actlist, nactive, i, j, bond, ma_type 1223 integer active(nw_max_atom) 1224c 1225 if (.not. rtdb_ma_get(rtdb, 'oniom:actlist', ma_type, 1226 $ nactive, l_actlist)) return 1227c 1228 if (.not. ma_get_index(l_actlist, k_actlist)) 1229 $ call errquit('grad_act_at: ma_get_inded failed',l_actlist, 1230 & MA_ERR) 1231c 1232c Prune out atoms not in the subset 1233c 1234 do i = 1, nat + nbonds 1235 active(i) = 0 1236 end do 1237 do i = 1, nactive 1238 j = int_mb(k_actlist+i-1) 1239 if (j.le.nat) active(j) = 1 1240 end do 1241c 1242c Ensure link atoms are active if necessary 1243c 1244 do bond = 1, nbonds 1245 i = bonds(1,bond) 1246 j = bonds(2,bond) 1247 if (active(i).gt.0 .or. active(j).gt.0) active(nat+bond) = 1 1248 end do 1249c 1250c Compress the active flags down into a list of atoms 1251c 1252 nactive = 0 1253 do i = 1, nat + nbonds 1254 if (active(i).gt.0) then 1255 nactive = nactive + 1 1256 active(nactive) = i 1257 end if 1258 end do 1259c 1260 if (nactive .eq. 0) call errquit('oniom: nactive=0?',0, 1261 & INPUT_ERR) 1262c 1263 if (.not. rtdb_put(rtdb, 'geometry:actlist', mt_int, nactive, 1264 $ active)) call errquit('oniom: failed storing actlist',0, 1265 & RTDB_ERR) 1266 if (.not. ma_free_heap(l_actlist)) call errquit 1267 $ ('oniom: free of actlist?',0, MA_ERR) 1268c 1269 end 1270 subroutine oniom_theory_input(theory,basname,ecpname,input) 1271 implicit none 1272#include "errquit.fh" 1273 character*(*) theory, basname, ecpname, input ! [output] 1274#include "inp.fh" 1275c 1276c (high||medium||low) theory [basis basname] [ecp ecpname] [input inpstr] 1277c 1278c Read oniom theory input ... defaults assumed set outside 1279c 1280c The theory string should always be converted to lower case as the 1281c task layer uses lower case strings only. 1282c 1283 character*128 test 1284c 1285 if (.not. inp_a(theory)) call errquit 1286 $ ('oniom_input: failed reading theory',0, INPUT_ERR) 1287 call inp_lcase(theory) 1288c 1289 10 if (.not. inp_a(test)) return 1290 if (inp_compare(.false.,'basis',test)) then 1291 if (.not. inp_a(basname)) call errquit 1292 $ ('oniom_input: failed reading basname',0, INPUT_ERR) 1293 else if (inp_compare(.false.,'ecp',test)) then 1294 if (.not. inp_a(ecpname)) call errquit 1295 $ ('oniom_input: failed reading ecpname',0, INPUT_ERR) 1296 else if (inp_compare(.false.,'input',test)) then 1297 if (.not. inp_a(input)) call errquit 1298 $ ('oniom_input: failed reading input [input]',0, INPUT_ERR) 1299 else 1300 call errquit('oniom: unknown keyword for high||medium||low',0, 1301 & INPUT_ERR) 1302 end if 1303 goto 10 1304c 1305 end 1306 subroutine oniom_param(rtdb, geom, basis, ecp, charge, input) 1307 implicit none 1308#include "errquit.fh" 1309 integer rtdb 1310 character*(*) geom, basis, ecp, input 1311 double precision charge 1312#include "rtdb.fh" 1313#include "mafdecls.fh" 1314c 1315 logical ignore 1316 logical nw_inp_from_character 1317 external nw_inp_from_character 1318c 1319 if (charge .ne. -999.0d0) then 1320 if (.not. rtdb_put(rtdb,'charge',mt_dbl,1,charge)) 1321 $ call errquit('oniom:setting charge?',0, RTDB_ERR) 1322 end if 1323 if (geom .ne. ' ') then 1324 if (.not. rtdb_cput(rtdb,'geometry',1,geom)) 1325 $ call errquit('oniom: setting geometry?',0, RTDB_ERR) 1326 end if 1327c 1328 if (basis .ne. ' ') then 1329 if (.not. rtdb_cput(rtdb, 'ao basis', 1, basis)) 1330 $ call errquit('oniom: setting basis?',0, RTDB_ERR) 1331 end if 1332c 1333 if (ecp .ne. ' ') then 1334 if (.not. rtdb_cput(rtdb, 'ecp basis', 1, ecp)) 1335 $ call errquit('oniom: setting ecp?',0, RTDB_ERR) 1336 else 1337 ignore = rtdb_delete(rtdb,'ecp basis') 1338 end if 1339c 1340 if (input .ne. ' ') then 1341 if (.not. nw_inp_from_character(rtdb,input)) 1342 $ call errquit('oniom: error processing input string',0, 1343 & INPUT_ERR) 1344 endif 1345c 1346 end 1347