1 subroutine qmmm_check_forces1(rtdb) 2 implicit none 3#include "errquit.fh" 4#include "global.fh" 5#include "mafdecls.fh" 6#include "stdio.fh" 7#include "util.fh" 8#include "rtdb.fh" 9#include "geom.fh" 10#include "qmmm.fh" 11#include "qmmm_params.fh" 12#include "mm_utils.fh" 13#include "inp.fh" 14c 15 integer rtdb 16c 17 integer ntot 18 integer i_c,h_c 19 integer i_c0,h_c0 20 integer i_ai,h_ai 21 integer i_g,h_g 22 double precision c0,gref 23 integer ia,i,k 24 double precision e0, ep, em, grad, hess, step 25 double precision step0 26 integer i1,k1 27c 28 character*32 pname 29 character*255 tag 30 character*30 region(3) 31 integer nregion 32c 33 pname = "qmmm_check_forces" 34c 35 if (ga_nodeid() .eq. 0) then 36 write(6,*) 37 write(6,*) '@ Checking forces' 38 write(6,*) 39 call util_flush(6) 40 end if 41c 42c 43c region definitions 44c ------------------ 45 tag ="qmmm:region" 46 if (.not.rtdb_get(rtdb,tag(1:inp_strlen(tag))//"_n", 47 > mt_int,1,nregion)) 48 > call errquit(pname//tag,0,RTDB_ERR) 49 if(nregion.gt.3) 50 > call errquit(pname//"too many regions",0,0) 51 if (.not.rtdb_cget(rtdb,tag,nregion,region)) 52 > call errquit(pname//tag,0,RTDB_ERR) 53c 54c define set of active atoms 55c -------------------------- 56 call qmmm_cons_reset() 57 call qmmm_cons_set("fix","solute") 58 call qmmm_cons_set("fix","solvent") 59 call qmmm_cons_set("free",region(1)) 60c 61 62c 63c get total number of atoms 64c ------------------------- 65c call qmmm_cons_get_nacts(ntot) 66 ntot = 1 67 if (ga_nodeid() .eq. 0) then 68 write(6,*) 69 write(6,*) '@ number of atoms',ntot 70 write(6,*) 71 call util_flush(6) 72 end if 73c 74 if(.not.ma_alloc_get(MT_INT, ntot, 'qmmm ref index array', 75 & h_ai,i_ai) ) call errquit( 76 & 'qmmm_data_alloc: unable to allocate heap space', 77 & ntot, MA_ERR) 78 call ifill(ntot,-1,int_mb(i_ai),1) 79 80 if(.not.ma_alloc_get(MT_DBL, 3*ntot, 'qmmm ref grad array', 81 & h_g, i_g) ) call errquit( 82 & 'qmmm_data_alloc: unable to allocate heap space', 83 & 3*ntot, MA_ERR) 84 call dfill(3*ntot,0,dbl_mb(i_g),1) 85 86 if(.not.ma_alloc_get(MT_DBL, 3*ntot, 'qmmm ref coord array', 87 & h_c,i_c) ) call errquit( 88 & 'qmmm_data_alloc: unable to allocate heap space', 89 & ntot, MA_ERR) 90 call dfill(3*ntot,0,dbl_mb(i_c),1) 91 92c call qmmm_cons_get_i_acts(int_mb(i_ai)) 93 int_mb(i_ai) = 10 94 call mm_get_solute_coord_raw(ntot, 95 > int_mb(i_ai), 96 > dbl_mb(i_c)) 97 98 call md_sp_qmmm() 99 call prp_print() 100 if (.not. rtdb_get(rtdb,'md:energy', mt_dbl, 1, e0)) 101 $ call errquit('driver: could not get energy',0, RTDB_ERR) 102 call mm_get_solute_force_raw(ntot, 103 & int_mb(i_ai), 104 & dbl_mb(i_g)) 105 106 if(.not.rtdb_get(rtdb,'qmmm:step',mt_dbl, 107 > 1,step0)) 108 > step0=0.01d0 109 110 call mm_print_system() 111 112c 113c loop over all atoms 114c ------------------- 115 if (ga_nodeid() .eq. 0) 116 > write(6,8) "comp","anal-g","num-g","error","step","ep","em" 117 8 format(1x,'checkgrad ',A5,4A16,2x,2A12) 118 119 i=0 120 do ia=1,ntot 121 do k=1,3 122 i=i+1 123 c0=dbl_mb(i_c+i-1) 124 gref=dbl_mb(i_g+i-1) 125 126 step = step0 12710 continue 128 dbl_mb(i_c+i-1)=c0+step 129 call mm_set_solute_coord_raw(ntot, 130 > int_mb(i_ai), 131 > dbl_mb(i_c)) 132 133 call md_sp_qmmm() 134 if (.not. rtdb_get(rtdb,'md:energy', mt_dbl, 1, ep)) 135 $ call errquit('driver: could not get energy',0, RTDB_ERR) 136 call mm_print_system() 137 call prp_print() 138 write(*,*) "c,c0,ep",dbl_mb(i_c+i-1),c0,ep 139 140c if (abs(ep-e0) .lt. 1e-6) then 141c write(6,*) ' Increasing the step ', ep,e0,ep-e0, step 142c step = step*10.0d0 143c goto 10 144c else if (abs(ep-e0) .gt. 1e-2) then 145c write(6,*) ' Decreasing the step ', ep,e0,ep-e0, step 146c step = step/3.0 147c if(step.le.0.00000001d0) goto 12 148c goto 10 149c end if 15012 continue 151 152 dbl_mb(i_c+i-1)=c0-step 153 call mm_set_solute_coord_raw(ntot, 154 > int_mb(i_ai), 155 > dbl_mb(i_c)) 156 157 call md_sp_qmmm() 158 if (.not. rtdb_get(rtdb,'md:energy', mt_dbl, 1, em)) 159 $ call errquit('driver: could not get energy',0, RTDB_ERR) 160 161 write(*,*) "c,c0,em",dbl_mb(i_c+i-1),c0,em 162 163 grad = (ep - em) / (2.0d0*step) 164c 165 if (ga_nodeid() .eq. 0) 166 > write(6,7) int_mb(i_ai+ia-1), 167 > gref, grad,abs(gref-grad), step, ep,em 168 7 format(1x,'checkgrad ',i5,2f16.8,e16.8,f16.8,2x,2f16.8) 169 dbl_mb(i_c+i-1)=c0 170 end do 171 end do 172 173 174 end 175 176 subroutine qmmm_check_forces(rtdb) 177 implicit none 178#include "errquit.fh" 179#include "global.fh" 180#include "mafdecls.fh" 181#include "stdio.fh" 182#include "util.fh" 183#include "rtdb.fh" 184#include "geom.fh" 185#include "qmmm.fh" 186#include "qmmm_params.fh" 187#include "mm_utils.fh" 188#include "inp.fh" 189c 190 integer rtdb 191c 192 integer ntot 193 integer i_c,h_c 194 integer i_c0,h_c0 195 integer i_ai,h_ai 196 integer i_g,h_g 197 double precision c0,gref 198 integer ia,i,k 199 double precision e0, ep, em, grad, hess, step 200 double precision step0 201 integer i1,k1 202c 203 logical qmmm_energy_gradient 204 logical task_qmmm_gradient 205 external qmmm_energy_gradient 206 external task_qmmm_gradient 207 logical task_qmmm_energy 208 external task_qmmm_energy 209 character*32 pname 210 character*255 tag 211 character*30 region(3) 212 integer nregion 213c 214 pname = "qmmm_check_forces" 215c 216 if (ga_nodeid() .eq. 0) then 217 write(6,*) 218 write(6,*) '@ Checking forces' 219 write(6,*) 220 call util_flush(6) 221 end if 222c 223c 224c region definitions 225c ------------------ 226 tag ="qmmm:region" 227 if (.not.rtdb_get(rtdb,tag(1:inp_strlen(tag))//"_n", 228 > mt_int,1,nregion)) 229 > call errquit(pname//tag,0,RTDB_ERR) 230 if(nregion.gt.3) 231 > call errquit(pname//"too many regions",0,0) 232 if (.not.rtdb_cget(rtdb,tag,nregion,region)) 233 > call errquit(pname//tag,0,RTDB_ERR) 234c 235c define set of active atoms 236c -------------------------- 237 call qmmm_cons_reset() 238 call qmmm_cons_set("fix","solute") 239 call qmmm_cons_set("fix","solvent") 240 call qmmm_cons_set("free",region(1)) 241c 242 243 call mm_print_system() 244c 245c get total number of atoms 246c ------------------------- 247 call qmmm_cons_get_nacts(ntot) 248 if (ga_nodeid() .eq. 0) then 249 write(6,*) 250 write(6,*) '@ number of atoms',ntot 251 write(6,*) 252 call util_flush(6) 253 end if 254c 255 if(.not.ma_alloc_get(MT_INT, ntot, 'qmmm ref index array', 256 & h_ai,i_ai) ) call errquit( 257 & 'qmmm_data_alloc: unable to allocate heap space', 258 & ntot, MA_ERR) 259 call ifill(ntot,-1,int_mb(i_ai),1) 260 261 if(.not.ma_alloc_get(MT_DBL, 3*ntot, 'qmmm ref grad array', 262 & h_g, i_g) ) call errquit( 263 & 'qmmm_data_alloc: unable to allocate heap space', 264 & 3*ntot, MA_ERR) 265 call dfill(3*ntot,0,dbl_mb(i_g),1) 266 267 if(.not.ma_alloc_get(MT_DBL, 3*ntot, 'qmmm ref coord array', 268 & h_c,i_c) ) call errquit( 269 & 'qmmm_data_alloc: unable to allocate heap space', 270 & ntot, MA_ERR) 271 call dfill(3*ntot,0,dbl_mb(i_c),1) 272 273 call qmmm_cons_get_i_acts(int_mb(i_ai)) 274 call mm_get_solute_coord(ntot, 275 > int_mb(i_ai), 276 > dbl_mb(i_c)) 277 278 call md_sp_qmmm() 279 if (.not. qmmm_energy_gradient(rtdb,.true.)) 280 $ call errquit(pname//'qmmm_energy_gradient failed' 281 $ ,0, GEOM_ERR) 282 call qmmm_energy_rtdb_push(rtdb) 283 if (.not. rtdb_get(rtdb,'qmmm:energy', mt_dbl, 1, e0)) 284 $ call errquit('driver: could not get energy',0, RTDB_ERR) 285 call qmmm_print_energy(rtdb) 286 call mm_get_solute_force(ntot, 287 & int_mb(i_ai), 288 & dbl_mb(i_g)) 289 290 if(.not.rtdb_get(rtdb,'qmmm:step',mt_dbl, 291 > 1,step0)) 292 > step0=0.01d0 293 294c 295c loop over all atoms 296c ------------------- 297 if (ga_nodeid() .eq. 0) 298 > write(6,8) "comp","anal-g","num-g","error","step","ep","em" 299 8 format(1x,'checkgrad ',A5,4A16,2x,2A12) 300 301 i=0 302 do ia=1,ntot 303 do k=1,3 304 i=i+1 305 c0=dbl_mb(i_c+i-1) 306 gref=dbl_mb(i_g+i-1) 307 308 step = step0 30910 continue 310 dbl_mb(i_c+i-1)=c0+step 311 call mm_set_solute_coord(ntot, 312 > int_mb(i_ai), 313 > dbl_mb(i_c)) 314 315 call md_sp_qmmm() 316 if (.not. qmmm_energy_gradient(rtdb,.false.)) 317 $ call errquit(pname//'qmmm_energy_gradient failed' 318 $ ,0, GEOM_ERR) 319 call qmmm_energy_rtdb_push(rtdb) 320 if (.not. rtdb_get(rtdb,'qmmm:energy', mt_dbl, 1, ep)) 321 $ call errquit('driver: could not get energy',0, RTDB_ERR) 322 call qmmm_print_energy(rtdb) 323 324c if (abs(ep-e0) .lt. 1e-6) then 325c write(6,*) ' Increasing the step ', ep,e0,ep-e0, step 326c step = step*10.0d0 327c goto 10 328c else if (abs(ep-e0) .gt. 1e-2) then 329c write(6,*) ' Decreasing the step ', ep,e0,ep-e0, step 330c step = step/3.0 331c if(step.le.0.00000001d0) goto 12 332c goto 10 333c end if 33412 continue 335 336 dbl_mb(i_c+i-1)=c0-step 337 call mm_set_solute_coord(ntot, 338 > int_mb(i_ai), 339 > dbl_mb(i_c)) 340 341 call md_sp_qmmm() 342 if (.not. qmmm_energy_gradient(rtdb,.false.)) 343 $ call errquit(pname//'qmmm_energy_gradient failed' 344 $ ,0, GEOM_ERR) 345 call qmmm_energy_rtdb_push(rtdb) 346 if (.not. rtdb_get(rtdb,'qmmm:energy', mt_dbl, 1, em)) 347 $ call errquit('driver: could not get energy',0, RTDB_ERR) 348 call qmmm_print_energy(rtdb) 349 350 351 grad = (ep - em) / (2.0d0*step) 352c 353 if (ga_nodeid() .eq. 0) 354 > write(6,7) int_mb(i_ai+ia-1), 355 > gref, grad,abs(gref-grad), step, ep,em 356 7 format(1x,'checkgrad ',i5,2f16.8,e16.8,f16.8,2x,2f16.8) 357 dbl_mb(i_c+i-1)=c0 358 end do 359 end do 360 361 362 end 363 364 subroutine qmmm_check_forces0(rtdb) 365 implicit none 366#include "errquit.fh" 367#include "global.fh" 368#include "mafdecls.fh" 369#include "stdio.fh" 370#include "util.fh" 371#include "rtdb.fh" 372#include "geom.fh" 373#include "qmmm.fh" 374#include "qmmm_params.fh" 375#include "mm_utils.fh" 376c 377 integer rtdb 378c 379 integer ntot 380 integer i_c,h_c 381 integer i_c0,h_c0 382 integer i_ai,h_ai 383 integer i_g,h_g 384 double precision c0,gref 385 integer ia,i,k 386 double precision e0, ep, em, grad, hess, step 387 double precision step0 388 integer i1,k1 389c 390 logical qmmm_energy_gradient 391 external qmmm_energy_gradient 392 logical task_qmmm_energy 393 external task_qmmm_energy 394 character*32 pname 395 396 pname = "qmmm_check_forces" 397c 398 call md_sp() 399 if (.not. rtdb_get(rtdb,'md:energy',mt_dbl,1,e0)) 400 $ call errquit('qmmm: failed put energy', 0, RTDB_ERR) 401 e0 = e0/cau2kj 402c 403 if (ga_nodeid() .eq. 0) then 404 write(6,*) 405 write(6,*) '@ Checking forces' 406 write(6,*) 407 call util_flush(6) 408 end if 409c 410c get total number of atoms 411c ------------------------- 412 call mm_get_solute_tot_na_gen(ntot,mm_link) 413 if (ga_nodeid() .eq. 0) then 414 write(6,*) 415 write(6,*) '@ number of atoms',ntot 416 write(6,*) 417 call util_flush(6) 418 end if 419 420c 421 if(.not.ma_alloc_get(MT_INT, ntot, 'qmmm ref index array', 422 & h_ai,i_ai) ) call errquit( 423 & 'qmmm_data_alloc: unable to allocate heap space', 424 & ntot, MA_ERR) 425 call ifill(ntot,-1,int_mb(i_ai),1) 426 427 if(.not.ma_alloc_get(MT_DBL, 3*ntot, 'qmmm ref grad array', 428 & h_g, i_g) ) call errquit( 429 & 'qmmm_data_alloc: unable to allocate heap space', 430 & 3*ntot, MA_ERR) 431 call dfill(3*ntot,0,dbl_mb(i_g),1) 432 433 if(.not.ma_alloc_get(MT_DBL, 3*ntot, 'qmmm ref coord array', 434 & h_c,i_c) ) call errquit( 435 & 'qmmm_data_alloc: unable to allocate heap space', 436 & ntot, MA_ERR) 437 call dfill(3*ntot,0,dbl_mb(i_c),1) 438 439 call mm_get_solute_ind_gen(ntot, 440 > mm_link, 441 > int_mb(i_ai)) 442 443 call mm_get_solute_coord_gen(ntot, 444 > mm_link, 445 > int_mb(i_ai), 446 > dbl_mb(i_c)) 447 448 call mm_get_solute_force_gen(ntot, 449 > mm_link, 450 > int_mb(i_ai), 451 > dbl_mb(i_g)) 452 453 if(.not.rtdb_get(rtdb,'qmmm:step',mt_dbl, 454 > 1,step0)) 455 > step0=0.01d0 456 457c do i=1,3*ntot 458c dbl_mb(i_g+i-1)=dbl_mb(i_g+i-1)*cau2kj/cau2nm 459c end do 460c 461c loop over all atoms 462c ------------------- 463 if (ga_nodeid() .eq. 0) 464 > write(6,8) "comp","anal-g","num-g","error","step","ep","em" 465 8 format(1x,'checkgrad ',A5,4A16,2x,2A12) 466 467 i=0 468 do ia=1,ntot 469 do k=1,3 470 i=i+1 471 c0=dbl_mb(i_c+i-1) 472 gref=dbl_mb(i_g+i-1) 473 474 step = step0 47510 continue 476 dbl_mb(i_c+i-1)=c0+step 477 call mm_set_solute_coord_gen(ntot, 478 > mm_link, 479 > int_mb(i_ai), 480 > dbl_mb(i_c)) 481 call md_sp() 482 if (.not. rtdb_get(rtdb,'md:energy',mt_dbl,1,ep)) 483 $ call errquit('qmmm: failed put energy', 0, RTDB_ERR) 484 ep = ep/cau2kj 485 486 if (abs(ep-e0) .lt. 1e-6) then 487 write(6,*) ' Increasing the step ', ep,e0,ep-e0, step 488 step = step*10.0d0 489 goto 10 490 else if (abs(ep-e0) .gt. 1e-2) then 491 write(6,*) ' Decreasing the step ', ep,e0,ep-e0, step 492 step = step/3.0 493 if(step.le.0.00000001d0) goto 12 494 goto 10 495 end if 49612 continue 497 498 dbl_mb(i_c+i-1)=c0-step 499 call mm_set_solute_coord_gen(ntot, 500 > mm_link, 501 > int_mb(i_ai), 502 > dbl_mb(i_c)) 503c 504 call md_sp() 505 if (.not. rtdb_get(rtdb,'md:energy',mt_dbl,1,em)) 506 $ call errquit('qmmm: failed put energy', 0, RTDB_ERR) 507 em = em/cau2kj 508c 509 grad = (ep - em) / (2.0d0*step) 510c 511 if (ga_nodeid() .eq. 0) 512 > write(6,7) int_mb(i_ai+ia-1), 513 > gref, grad,abs(gref-grad), step, ep,em 514 7 format(1x,'checkgrad ',i5,2f16.8,e16.8,f16.8,2x,2f16.8) 515 dbl_mb(i_c+i-1)=c0 516 end do 517 end do 518 519 520 end 521 522c $Id$ 523