1c $I#d: cons.F,v 1.1 2004/01/28 01:30:59 marat Exp $ 2 subroutine cons_unload_hbonds() 3 implicit none 4#include "errquit.fh" 5#include "mafdecls.fh" 6#include "cons.fh" 7 integer nhb,h_rhb,h_khb,h_ijhb 8 9 h_ijhb = cons_get_h_hbond_id() 10 h_khb = cons_get_h_hbond_k() 11 h_rhb = cons_get_h_hbond_r() 12c 13c unload harmonic constraints if any 14c 15 call cons_get_hbond_nhb(nhb) 16 if(nhb.gt.0) then 17 if (.not. ma_free_heap(h_rhb) ) call errquit( 18 & 'cons_bond_input: unable to free h_rhb', 19 & 0, MA_ERR) 20 if (.not. ma_free_heap(h_khb) ) call errquit( 21 & 'cons_bond_input: unable to free h_khb', 22 & 0, MA_ERR) 23 if (.not. ma_free_heap(h_ijhb) ) call errquit( 24 & 'cons_bond_input: unable to free h_ijhb', 25 & 0, MA_ERR) 26 call cons_set_hbond_nhb(0) 27 end if 28 end 29c 30 subroutine cons_load_hbonds(namespace,rtdb) 31 implicit none 32#include "errquit.fh" 33#include "mafdecls.fh" 34#include "rtdb.fh" 35#include "util.fh" 36#include "geom.fh" 37c 38 character*(*) namespace 39 integer rtdb 40c 41 logical status 42 integer nhb 43 integer i_rhb,i_khb,i_ijhb 44 integer h_rhb,h_khb,h_ijhb 45 character*255 tag_id 46 character*255 tag_n 47 character*255 tag_r 48 character*255 tag_k 49 integer ma_type,ma_n 50c 51 call cons_unload_hbonds() 52c 53 call cons_hbond_id_tag(namespace,tag_id) 54 call cons_hbond_n_tag(namespace,tag_n) 55 call cons_hbond_k_tag(namespace,tag_k) 56 call cons_hbond_r_tag(namespace,tag_r) 57c 58c load harmonic constraints 59c 60 status=rtdb_get(rtdb,tag_n,mt_int,1,nhb) 61c 62c exit if no harm bonds 63 if(.not.status .or. nhb.eq.0) then 64c call errquit( 65c > 'cons_load_hbonds: unable to get number of harm bonds:', 66c > nhb, MA_ERR) 67 return 68 end if 69 70 if ( .not. rtdb_ma_get( rtdb, tag_id,ma_type, ma_n, 71 & h_ijhb) ) call errquit( 72 & 'cons_load_hbonds: unable to allocate cons scratch space', 73 & 2*nhb, MA_ERR) 74 if ( .not. ma_get_index(h_ijhb, 75 & i_ijhb) ) call errquit( 76 & 'cons_load_hbonds: unable to allocate cons scratch space', 77 & 2*nhb, MA_ERR) 78 79 if ( .not. rtdb_ma_get( rtdb, tag_k,ma_type, ma_n, 80 & h_khb) ) call errquit( 81 & 'cons_load_hbonds: unable to allocate cons scratch space', 82 & 2*nhb, MA_ERR) 83 if ( .not. ma_get_index(h_khb, 84 & i_khb) ) call errquit( 85 & 'cons_load_hbonds: unable to allocate cons scratch space', 86 & 2*nhb, MA_ERR) 87 88 if ( .not. rtdb_ma_get( rtdb, tag_r,ma_type, ma_n, 89 & h_rhb) ) call errquit( 90 & 'cons_load_hbonds: unable to allocate cons scratch space', 91 & 2*nhb, MA_ERR) 92 if ( .not. ma_get_index(h_rhb, 93 & i_rhb) ) call errquit( 94 & 'cons_load_hbonds: unable to allocate cons scratch space', 95 & 2*nhb, MA_ERR) 96 97 98 call cons_set_hbond_nhb(nhb) 99 call cons_set_h_hbond_id(h_ijhb) 100 call cons_set_h_hbond_k(h_khb) 101 call cons_set_h_hbond_r(h_rhb) 102 103 end 104c 105 subroutine cons_delete_hbonds(namespace,rtdb) 106 implicit none 107#include "errquit.fh" 108#include "rtdb.fh" 109c 110 character*(*) namespace 111 integer rtdb 112c 113 logical status 114 character*255 tag_id 115 character*255 tag_n 116 character*255 tag_r 117 character*255 tag_k 118 119 call cons_hbond_id_tag(namespace,tag_id) 120 call cons_hbond_n_tag(namespace,tag_n) 121 call cons_hbond_k_tag(namespace,tag_k) 122 call cons_hbond_r_tag(namespace,tag_r) 123 124 status = rtdb_delete(rtdb,tag_n) 125 status = status .and. rtdb_delete(rtdb,tag_id) 126 status = status .and. rtdb_delete(rtdb,tag_k) 127 status = status .and. rtdb_delete(rtdb,tag_r) 128 129 end 130c 131 subroutine cons_add_hbond_energy(rtdb,energy) 132 implicit none 133#include "errquit.fh" 134#include "global.fh" 135#include "mafdecls.fh" 136#include "rtdb.fh" 137#include "util.fh" 138#include "stdio.fh" 139#include "cons.fh" 140c 141 integer rtdb 142 double precision energy 143 logical status 144 logical oprint 145 integer inb 146 integer iat,jat 147 double precision r 148 double precision r0,k 149 double precision e,enb 150 integer i_c,i_rhb,i_khb,i_ijhb 151 integer nhb 152c 153 call util_print_push() 154 call util_print_rtdb_load(rtdb,'cons') 155 oprint = util_print('energy', print_default) 156 oprint =oprint .and. (ga_nodeid().eq.0) 157c 158 call cons_get_hbond_nhb(nhb) 159 i_c = cons_get_i_c() 160 i_ijhb = cons_get_i_hbond_id() 161 i_khb = cons_get_i_hbond_k() 162 i_rhb = cons_get_i_hbond_r() 163 164 e=0.0d0 165 do inb=1,nhb 166 iat=int_mb(i_ijhb+2*(inb-1)) 167 jat=int_mb(i_ijhb+2*(inb-1)+1) 168 r0 =dbl_mb(i_rhb+inb-1) 169 k =dbl_mb(i_khb+inb-1) 170 171 call cons_spring_energy(k,r0, 172 > dbl_mb(i_c+(iat-1)*3), 173 > dbl_mb(i_c+(jat-1)*3), 174 > r, 175 > enb) 176 177 if(oprint) then 178 write(luout,*)"adding spring # ",inb 179 write(luout,*)" spring parameters (i,j,k,r0):",iat,jat,k,r0 180 write(luout,*)" spring length and energy :",r,enb 181 end if 182 e=e+enb 183 end do 184 if(rtdb_get(rtdb, 'cons:simulate', mt_log, 1, status)) then 185 write(luout,*) "cons energy simulation" 186 energy=e 187 else 188 energy = energy + e 189 end if 190c 191 call util_print_pop() 192 return 193 end 194c 195 196 subroutine cons_unload_hbondings() 197 implicit none 198#include "errquit.fh" 199#include "mafdecls.fh" 200#include "cons.fh" 201 integer nhc 202 integer h_hbondings_n0 203 integer h_hbondings_indx 204 integer h_hbondings_coef 205 integer h_hbondings_k0 206 integer h_hbondings_gamma0 207 208 h_hbondings_n0 = cons_get_h_hbondings_n0() 209 h_hbondings_indx = cons_get_h_hbondings_indx() 210 h_hbondings_coef = cons_get_h_hbondings_coef() 211 h_hbondings_k0 = cons_get_h_hbondings_k0() 212 h_hbondings_gamma0 = cons_get_h_hbondings_gamma0() 213c 214c unload harmonic constraints if any 215c 216 call cons_get_hbondings_nhc(nhc) 217 if(nhc.gt.0) then 218 if (.not. ma_free_heap(h_hbondings_n0) ) call errquit( 219 & 'cons_bonding_input: unable to free h_hbondings_n0', 220 & 0, MA_ERR) 221 if (.not. ma_free_heap(h_hbondings_indx) ) call errquit( 222 & 'cons_bonding_input: unable to free h_hbondings_indx', 223 & 0, MA_ERR) 224 if (.not. ma_free_heap(h_hbondings_coef) ) call errquit( 225 & 'cons_bonding_input: unable to free h_hbondings_coef', 226 & 0, MA_ERR) 227 if (.not. ma_free_heap(h_hbondings_k0) ) call errquit( 228 & 'cons_bonding_input: unable to free h_hbondings_k0', 229 & 0, MA_ERR) 230 if (.not. ma_free_heap(h_hbondings_gamma0) ) call errquit( 231 & 'cons_bonding_input: unable to free h_hbondings_gamma0', 232 & 0, MA_ERR) 233 call cons_set_hbondings_nhc(0) 234 end if 235 end 236 237c 238 subroutine cons_load_hbondings(namespace,rtdb) 239 implicit none 240#include "errquit.fh" 241#include "mafdecls.fh" 242#include "rtdb.fh" 243#include "util.fh" 244#include "geom.fh" 245c 246 character*(*) namespace 247 integer rtdb 248c 249 logical status 250 integer nhc 251 integer h_hbondings_n0,i_hbondings_n0 252 integer h_hbondings_indx,i_hbondings_indx 253 integer h_hbondings_coef,i_hbondings_coef 254 integer h_hbondings_k0,i_hbondings_k0 255 integer h_hbondings_gamma0,i_hbondings_gamma0 256 character*255 tag_n 257 character*255 tag_n0 258 character*255 tag_indx 259 character*255 tag_coef 260 character*255 tag_k0 261 character*255 tag_gamma0 262 integer ma_type,ma_n 263c 264 call cons_unload_hbondings() 265c 266 call cons_hbondings_n_tag(namespace,tag_n) 267 268 call cons_hbondings_n0_tag(namespace,tag_n0) 269 call cons_hbondings_indx_tag(namespace,tag_indx) 270 call cons_hbondings_coef_tag(namespace,tag_coef) 271 call cons_hbondings_k0_tag(namespace,tag_k0) 272 call cons_hbondings_gamma0_tag(namespace,tag_gamma0) 273c 274c load harmonic constraints 275c 276 status=rtdb_get(rtdb,tag_n,mt_int,1,nhc) 277c 278c exit if no harm bonds 279 if(.not.status .or. nhc.eq.0) then 280c call errquit( 281c > 'cons_load_hbonds: unable to get number of harm bonds:', 282c > nhb, MA_ERR) 283 return 284 end if 285 286 if ( .not. rtdb_ma_get( rtdb, tag_n0,ma_type, ma_n, 287 & h_hbondings_n0) ) call errquit( 288 & 'cons_load_hbonds: unable to allocate cons scratch space', 289 & 2*nhc, MA_ERR) 290 if ( .not. ma_get_index(h_hbondings_n0, 291 & i_hbondings_n0) ) call errquit( 292 & 'cons_load_hbonds: unable to allocate cons scratch space', 293 & 2*nhc, MA_ERR) 294 295 if ( .not. rtdb_ma_get( rtdb, tag_indx,ma_type, ma_n, 296 & h_hbondings_indx) ) call errquit( 297 & 'cons_load_hbonds: unable to allocate cons scratch space', 298 & 2*nhc, MA_ERR) 299 if ( .not. ma_get_index(h_hbondings_indx, 300 & i_hbondings_indx) ) call errquit( 301 & 'cons_load_hbonds: unable to allocate cons scratch space', 302 & 2*nhc, MA_ERR) 303 304 if ( .not. rtdb_ma_get( rtdb, tag_coef,ma_type, ma_n, 305 & h_hbondings_coef) ) call errquit( 306 & 'cons_load_hbonds: unable to allocate cons scratch space', 307 & 2*nhc, MA_ERR) 308 if ( .not. ma_get_index(h_hbondings_coef, 309 & i_hbondings_coef) ) call errquit( 310 & 'cons_load_hbonds: unable to allocate cons scratch space', 311 & 2*nhc, MA_ERR) 312 313 if ( .not. rtdb_ma_get( rtdb, tag_k0,ma_type, ma_n, 314 & h_hbondings_k0) ) call errquit( 315 & 'cons_load_hbonds: unable to allocate cons scratch space', 316 & 2*nhc, MA_ERR) 317 if ( .not. ma_get_index(h_hbondings_k0, 318 & i_hbondings_k0) ) call errquit( 319 & 'cons_load_hbonds: unable to allocate cons scratch space', 320 & 2*nhc, MA_ERR) 321 322 if ( .not. rtdb_ma_get( rtdb, tag_gamma0,ma_type, ma_n, 323 & h_hbondings_gamma0) ) call errquit( 324 & 'cons_load_hbonds: unable to allocate cons scratch space', 325 & 2*nhc, MA_ERR) 326 if ( .not. ma_get_index(h_hbondings_gamma0, 327 & i_hbondings_gamma0) ) call errquit( 328 & 'cons_load_hbonds: unable to allocate cons scratch space', 329 & 2*nhc, MA_ERR) 330 331 332 call cons_set_hbondings_nhc(nhc) 333 call cons_set_hbondings_n0(h_hbondings_n0) 334 call cons_set_hbondings_indx(h_hbondings_indx) 335 call cons_set_hbondings_coef(h_hbondings_coef) 336 call cons_set_hbondings_k0(h_hbondings_k0) 337 call cons_set_hbondings_gamma0(h_hbondings_gamma0) 338 339 return 340 end 341 342c 343 subroutine cons_delete_hbondings(namespace,rtdb) 344 implicit none 345#include "errquit.fh" 346#include "rtdb.fh" 347c 348 character*(*) namespace 349 integer rtdb 350c 351 logical status 352 character*255 tag_n 353 character*255 tag_n0 354 character*255 tag_indx 355 character*255 tag_coef 356 character*255 tag_k0 357 character*255 tag_gamma0 358 359 call cons_hbondings_n_tag(namespace,tag_n) 360 call cons_hbondings_n0_tag(namespace,tag_n0) 361 call cons_hbondings_indx_tag(namespace,tag_indx) 362 call cons_hbondings_coef_tag(namespace,tag_coef) 363 call cons_hbondings_k0_tag(namespace,tag_k0) 364 call cons_hbondings_k0_tag(namespace,tag_gamma0) 365 366 status = rtdb_delete(rtdb,tag_n) 367 status = status .and. rtdb_delete(rtdb,tag_n0) 368 status = status .and. rtdb_delete(rtdb,tag_indx) 369 status = status .and. rtdb_delete(rtdb,tag_coef) 370 status = status .and. rtdb_delete(rtdb,tag_k0) 371 status = status .and. rtdb_delete(rtdb,tag_gamma0) 372 373 end 374 375 376c 377 subroutine cons_add_hbondings_energy(rtdb,energy) 378 implicit none 379#include "errquit.fh" 380#include "global.fh" 381#include "mafdecls.fh" 382#include "rtdb.fh" 383#include "util.fh" 384#include "stdio.fh" 385#include "cons.fh" 386c 387 integer rtdb 388 double precision energy 389 logical oprint 390 integer i,ii,n,j 391 double precision K0,gamma0,gamm 392 double precision e,enb 393 integer i_c,i_n,i_indx,i_coef,i_k0,i_gamma0 394 integer nhc 395 396c 397 call util_print_push() 398 call util_print_rtdb_load(rtdb,'cons') 399 oprint = util_print('energy', print_default) 400 oprint =oprint .and. (ga_nodeid().eq.0) 401c 402 call cons_get_hbondings_nhc(nhc) 403 i_c = cons_get_i_c() 404 i_n = cons_get_i_hbondings_n0() 405 i_indx = cons_get_i_hbondings_indx() 406 i_coef = cons_get_i_hbondings_coef() 407 i_k0 = cons_get_i_hbondings_k0() 408 i_gamma0 = cons_get_i_hbondings_gamma0() 409 410 e=0.0d0 411 ii = 0 412 do i=1,nhc 413 n = int_mb(i_n+i-1) 414 K0 = dbl_mb(i_k0+i-1) 415 gamma0 = dbl_mb(i_gamma0+i-1) 416 417 call cons_bondings_energy(n,int_mb(i_indx+2*ii), 418 > dbl_mb(i_coef+ii), 419 > K0,gamma0,dbl_mb(i_c),gamm,enb) 420 if(oprint) then 421 write(luout,*) 422 write(luout,'(A,I4)') "Adding bondings spring # ",i 423 write(luout,'(A,F18.9,E14.6)') 424 > " - spring parameters (K0,gamma0):", 425 > K0,gamma0 426 write(luout,'(A,E18.6)') 427 > " - gamma :",gamm 428 write(luout,'(A,F18.9)') 429 > " - spring energy :",enb 430 write(luout,'(A,F18.9)') 431 > " - previous/base energy :",energy 432 write(luout,'(A,F18.9)') 433 > " - constrained energy :",energy+enb 434 write(luout,*) " - coefficient index1 index2" 435 do j=1,n 436 write(luout,'(F16.6,I7,I7)') dbl_mb(i_coef+ii+j-1), 437 > int_mb(i_indx+2*ii+2*(j-1)), 438 > int_mb(i_indx+2*ii+2*(j-1)+1) 439 end do 440 write(luout,*) 441 end if 442 ii = ii + n 443 e=e+enb 444 end do 445 energy = energy + e 446c 447 call util_print_pop() 448 return 449 end 450 451 452c 453 subroutine cons_add_hbondings_egrad(rtdb,energy,gx) 454 implicit none 455#include "errquit.fh" 456#include "mafdecls.fh" 457#include "global.fh" 458#include "rtdb.fh" 459#include "util.fh" 460#include "stdio.fh" 461#include "cons.fh" 462c 463 integer rtdb 464 double precision gx(*) 465 double precision energy 466c local variables 467c --------------- 468 logical oprint 469 integer i,j,n,ii 470 integer nhc 471 double precision gamma0,K0,gamm 472 double precision e,enb 473 integer i_c,i_n,i_indx,i_coef,i_k0,i_gamma0 474c 475 call util_print_push() 476 call util_print_rtdb_load(rtdb,'cons') 477 oprint = util_print('energy', print_default) 478 oprint =oprint .and. (ga_nodeid().eq.0) 479c 480 call cons_get_hbondings_nhc(nhc) 481 i_c = cons_get_i_c() 482 i_n = cons_get_i_hbondings_n0() 483 i_indx = cons_get_i_hbondings_indx() 484 i_coef = cons_get_i_hbondings_coef() 485 i_k0 = cons_get_i_hbondings_k0() 486 i_gamma0 = cons_get_i_hbondings_gamma0() 487 488 489 e=0.0d0 490 ii = 0 491 do i=1,nhc 492 n = int_mb(i_n+i-1) 493 K0 = dbl_mb(i_k0+i-1) 494 gamma0 = dbl_mb(i_gamma0+i-1) 495 496 call cons_bondings_force(n,int_mb(i_indx+2*ii), 497 > dbl_mb(i_coef+ii), 498 > K0,gamma0,dbl_mb(i_c),gamm,enb,gx) 499 if(oprint) then 500 write(luout,*) 501 write(luout,'(A,I4)') "Adding bondings spring # ",i 502 write(luout,'(A,F18.9,E14.6)') 503 > " - spring parameters (K0,gamma0):", 504 > K0,gamma0 505 write(luout,'(A,E18.6)') 506 > " - gamma :",gamm 507 write(luout,'(A,F18.9)') 508 > " - spring energy :",enb 509 write(luout,'(A,F18.9)') 510 > " - previous/base energy :",energy 511 write(luout,'(A,F18.9)') 512 > " - constrained energy :",energy+enb 513 write(luout,*) " - coefficient index1 index2" 514 do j=1,n 515 write(luout,'(F16.6,I7,I7)') dbl_mb(i_coef+ii+j-1), 516 > int_mb(i_indx+2*ii+2*(j-1)), 517 > int_mb(i_indx+2*ii+2*(j-1)+1) 518 end do 519 write(luout,*) 520 end if 521 522 ii = ii + n 523 e=e+enb 524 end do 525 energy = energy + e 526c 527 call util_print_pop() 528 return 529 end 530 531 532 533 subroutine cons_unload_pbondings() 534 implicit none 535#include "errquit.fh" 536#include "mafdecls.fh" 537#include "cons.fh" 538 integer nhp 539 integer h_pbondings_n0 540 integer h_pbondings_indx 541 integer h_pbondings_coef 542 integer h_pbondings_k0 543 integer h_pbondings_gcut0 544 integer h_pbondings_gamma0 545 546 h_pbondings_n0 = cons_get_h_pbondings_n0() 547 h_pbondings_indx = cons_get_h_pbondings_indx() 548 h_pbondings_coef = cons_get_h_pbondings_coef() 549 h_pbondings_k0 = cons_get_h_pbondings_k0() 550 h_pbondings_gcut0 = cons_get_h_pbondings_gcut0() 551 h_pbondings_gamma0 = cons_get_h_pbondings_gamma0() 552c 553c unload harmonic constraints if any 554c 555 call cons_get_pbondings_nhp(nhp) 556 if(nhp.gt.0) then 557 if (.not. ma_free_heap(h_pbondings_n0) ) call errquit( 558 & 'cons_bonding_input: unable to free h_pbondings_n0', 559 & 0, MA_ERR) 560 if (.not. ma_free_heap(h_pbondings_indx) ) call errquit( 561 & 'cons_bonding_input: unable to free h_pbondings_indx', 562 & 0, MA_ERR) 563 if (.not. ma_free_heap(h_pbondings_coef) ) call errquit( 564 & 'cons_bonding_input: unable to free h_pbondings_coef', 565 & 0, MA_ERR) 566 if (.not. ma_free_heap(h_pbondings_k0) ) call errquit( 567 & 'cons_bonding_input: unable to free h_pbondings_k0', 568 & 0, MA_ERR) 569 if (.not. ma_free_heap(h_pbondings_gcut0) ) call errquit( 570 & 'cons_bonding_input: unable to free h_pbondings_gcut0', 571 & 0, MA_ERR) 572 if (.not. ma_free_heap(h_pbondings_gamma0) ) call errquit( 573 & 'cons_bonding_input: unable to free h_pbondings_gamma0', 574 & 0, MA_ERR) 575 call cons_set_pbondings_nhp(0) 576 end if 577 end 578 579c 580 subroutine cons_load_pbondings(namespace,rtdb) 581 implicit none 582#include "errquit.fh" 583#include "mafdecls.fh" 584#include "rtdb.fh" 585#include "util.fh" 586#include "geom.fh" 587c 588 character*(*) namespace 589 integer rtdb 590c 591 logical status 592 integer nhp 593 integer h_pbondings_n0,i_pbondings_n0 594 integer h_pbondings_indx,i_pbondings_indx 595 integer h_pbondings_coef,i_pbondings_coef 596 integer h_pbondings_k0,i_pbondings_k0 597 integer h_pbondings_gcut0,i_pbondings_gcut0 598 integer h_pbondings_gamma0,i_pbondings_gamma0 599 character*255 tag_n 600 character*255 tag_n0 601 character*255 tag_indx 602 character*255 tag_coef 603 character*255 tag_k0 604 character*255 tag_gcut0 605 character*255 tag_gamma0 606 integer ma_type,ma_n 607c 608 call cons_unload_pbondings() 609c 610 call cons_pbondings_n_tag(namespace,tag_n) 611 612 call cons_pbondings_n0_tag(namespace,tag_n0) 613 call cons_pbondings_indx_tag(namespace,tag_indx) 614 call cons_pbondings_coef_tag(namespace,tag_coef) 615 call cons_pbondings_k0_tag(namespace,tag_k0) 616 call cons_pbondings_gcut0_tag(namespace,tag_gcut0) 617 call cons_pbondings_gamma0_tag(namespace,tag_gamma0) 618c 619c load harmonic constraints 620c 621 status=rtdb_get(rtdb,tag_n,mt_int,1,nhp) 622c 623c exit if no harm bonds 624 if(.not.status .or. nhp.eq.0) then 625c call errquit( 626c > 'cons_load_hbonds: unable to get number of harm bonds:', 627c > nhb, MA_ERR) 628 return 629 end if 630 631 if ( .not. rtdb_ma_get( rtdb, tag_n0,ma_type, ma_n, 632 & h_pbondings_n0) ) call errquit( 633 & 'cons_load_hbonds: unable to allocate cons scratch space', 634 & 2*nhp, MA_ERR) 635 if ( .not. ma_get_index(h_pbondings_n0, 636 & i_pbondings_n0) ) call errquit( 637 & 'cons_load_hbonds: unable to allocate cons scratch space', 638 & 2*nhp, MA_ERR) 639 640 if ( .not. rtdb_ma_get( rtdb, tag_indx,ma_type, ma_n, 641 & h_pbondings_indx) ) call errquit( 642 & 'cons_load_hbonds: unable to allocate cons scratch space', 643 & 2*nhp, MA_ERR) 644 if ( .not. ma_get_index(h_pbondings_indx, 645 & i_pbondings_indx) ) call errquit( 646 & 'cons_load_hbonds: unable to allocate cons scratch space', 647 & 2*nhp, MA_ERR) 648 649 if ( .not. rtdb_ma_get( rtdb, tag_coef,ma_type, ma_n, 650 & h_pbondings_coef) ) call errquit( 651 & 'cons_load_hbonds: unable to allocate cons scratch space', 652 & 2*nhp, MA_ERR) 653 if ( .not. ma_get_index(h_pbondings_coef, 654 & i_pbondings_coef) ) call errquit( 655 & 'cons_load_hbonds: unable to allocate cons scratch space', 656 & 2*nhp, MA_ERR) 657 658 if ( .not. rtdb_ma_get( rtdb, tag_k0,ma_type, ma_n, 659 & h_pbondings_k0) ) call errquit( 660 & 'cons_load_hbonds: unable to allocate cons scratch space', 661 & 2*nhp, MA_ERR) 662 if ( .not. ma_get_index(h_pbondings_k0, 663 & i_pbondings_k0) ) call errquit( 664 & 'cons_load_hbonds: unable to allocate cons scratch space', 665 & 2*nhp, MA_ERR) 666 667 if ( .not. rtdb_ma_get( rtdb, tag_gcut0,ma_type, ma_n, 668 & h_pbondings_gcut0) ) call errquit( 669 & 'cons_load_hbonds: unable to allocate cons scratch space', 670 & 2*nhp, MA_ERR) 671 if ( .not. ma_get_index(h_pbondings_gcut0, 672 & i_pbondings_gcut0) ) call errquit( 673 & 'cons_load_hbonds: unable to allocate cons scratch space', 674 & 2*nhp, MA_ERR) 675 676 if ( .not. rtdb_ma_get( rtdb, tag_gamma0,ma_type, ma_n, 677 & h_pbondings_gamma0) ) call errquit( 678 & 'cons_load_hbonds: unable to allocate cons scratch space', 679 & 2*nhp, MA_ERR) 680 if ( .not. ma_get_index(h_pbondings_gamma0, 681 & i_pbondings_gamma0) ) call errquit( 682 & 'cons_load_hbonds: unable to allocate cons scratch space', 683 & 2*nhp, MA_ERR) 684 685 686 call cons_set_pbondings_nhp(nhp) 687 call cons_set_pbondings_n0(h_pbondings_n0) 688 call cons_set_pbondings_indx(h_pbondings_indx) 689 call cons_set_pbondings_coef(h_pbondings_coef) 690 call cons_set_pbondings_k0(h_pbondings_k0) 691 call cons_set_pbondings_gcut0(h_pbondings_gcut0) 692 call cons_set_pbondings_gamma0(h_pbondings_gamma0) 693 694 return 695 end 696 697c 698 subroutine cons_delete_pbondings(namespace,rtdb) 699 implicit none 700#include "errquit.fh" 701#include "rtdb.fh" 702c 703 character*(*) namespace 704 integer rtdb 705c 706 logical status 707 character*255 tag_n 708 character*255 tag_n0 709 character*255 tag_indx 710 character*255 tag_coef 711 character*255 tag_k0 712 character*255 tag_gcut0 713 character*255 tag_gamma0 714 715 call cons_pbondings_n_tag(namespace,tag_n) 716 call cons_pbondings_n0_tag(namespace,tag_n0) 717 call cons_pbondings_indx_tag(namespace,tag_indx) 718 call cons_pbondings_coef_tag(namespace,tag_coef) 719 call cons_pbondings_k0_tag(namespace,tag_k0) 720 call cons_pbondings_gcut0_tag(namespace,tag_gcut0) 721 call cons_pbondings_gamma0_tag(namespace,tag_gamma0) 722 723 status = rtdb_delete(rtdb,tag_n) 724 status = status .and. rtdb_delete(rtdb,tag_n0) 725 status = status .and. rtdb_delete(rtdb,tag_indx) 726 status = status .and. rtdb_delete(rtdb,tag_coef) 727 status = status .and. rtdb_delete(rtdb,tag_k0) 728 status = status .and. rtdb_delete(rtdb,tag_gcut0) 729 status = status .and. rtdb_delete(rtdb,tag_gamma0) 730 731 end 732 733 734c 735 subroutine cons_add_pbondings_energy(rtdb,energy) 736 implicit none 737#include "errquit.fh" 738#include "global.fh" 739#include "mafdecls.fh" 740#include "rtdb.fh" 741#include "util.fh" 742#include "stdio.fh" 743#include "cons.fh" 744c 745 integer rtdb 746 double precision energy 747 logical oprint 748 integer i,ii,n,j 749 double precision K0,gcut0,gamma0,gamm 750 double precision e,enb 751 integer i_c,i_n,i_indx,i_coef,i_k0,i_gcut0,i_gamma0 752 integer nhp 753 754c 755 call util_print_push() 756 call util_print_rtdb_load(rtdb,'cons') 757 oprint = util_print('energy', print_default) 758 oprint =oprint .and. (ga_nodeid().eq.0) 759c 760 call cons_get_pbondings_nhp(nhp) 761 i_c = cons_get_i_c() 762 i_n = cons_get_i_pbondings_n0() 763 i_indx = cons_get_i_pbondings_indx() 764 i_coef = cons_get_i_pbondings_coef() 765 i_k0 = cons_get_i_pbondings_k0() 766 i_gcut0 = cons_get_i_pbondings_gcut0() 767 i_gamma0 = cons_get_i_pbondings_gamma0() 768 769 e=0.0d0 770 ii = 0 771 do i=1,nhp 772 n = int_mb(i_n+i-1) 773 K0 = dbl_mb(i_k0+i-1) 774 gcut0 = dbl_mb(i_gcut0+i-1) 775 gamma0 = dbl_mb(i_gamma0+i-1) 776 777 call cons_pbondings_energy(n,int_mb(i_indx+2*ii), 778 > dbl_mb(i_coef+ii), 779 > K0,gcut0,gamma0,dbl_mb(i_c),gamm,enb) 780 if(oprint) then 781 write(luout,*) 782 write(luout,'(A,I4)') "Adding bondings spring # ",i 783 write(luout,'(A,F18.9,F18.9,E14.6)') 784 > " - spring parameters (K0,gcut0,gamma0):", 785 > K0,gcut0,gamma0 786 write(luout,'(A,E18.6)') 787 > " - gamma :",gamm 788 write(luout,'(A,F18.9)') 789 > " - spring energy :",enb 790 write(luout,'(A,F18.9)') 791 > " - previous/base energy :",energy 792 write(luout,'(A,F18.9)') 793 > " - constrained energy :",energy+enb 794 write(luout,*) " - coefficient index1 index2" 795 do j=1,n 796 write(luout,'(F16.6,I7,I7)') dbl_mb(i_coef+ii+j-1), 797 > int_mb(i_indx+2*ii+2*(j-1)), 798 > int_mb(i_indx+2*ii+2*(j-1)+1) 799 end do 800 write(luout,*) 801 end if 802 ii = ii + n 803 e=e+enb 804 end do 805 energy = energy + e 806c 807 call util_print_pop() 808 return 809 end 810 811 812c 813 subroutine cons_add_pbondings_egrad(rtdb,energy,gx) 814 implicit none 815#include "errquit.fh" 816#include "mafdecls.fh" 817#include "global.fh" 818#include "rtdb.fh" 819#include "util.fh" 820#include "stdio.fh" 821#include "cons.fh" 822c 823 integer rtdb 824 double precision gx(*) 825 double precision energy 826c local variables 827c --------------- 828 logical oprint 829 integer i,j,n,ii 830 integer nhp 831 double precision gamma0,K0,gcut0,gamm 832 double precision e,enb 833 integer i_c,i_n,i_indx,i_coef,i_k0,i_gcut0,i_gamma0 834c 835 call util_print_push() 836 call util_print_rtdb_load(rtdb,'cons') 837 oprint = util_print('energy', print_default) 838 oprint =oprint .and. (ga_nodeid().eq.0) 839c 840 call cons_get_pbondings_nhp(nhp) 841 i_c = cons_get_i_c() 842 i_n = cons_get_i_pbondings_n0() 843 i_indx = cons_get_i_pbondings_indx() 844 i_coef = cons_get_i_pbondings_coef() 845 i_k0 = cons_get_i_pbondings_k0() 846 i_gcut0 = cons_get_i_pbondings_gcut0() 847 i_gamma0 = cons_get_i_pbondings_gamma0() 848 849 850 e=0.0d0 851 ii = 0 852 do i=1,nhp 853 n = int_mb(i_n+i-1) 854 K0 = dbl_mb(i_k0+i-1) 855 gcut0 = dbl_mb(i_gcut0+i-1) 856 gamma0 = dbl_mb(i_gamma0+i-1) 857 858 call cons_pbondings_force(n,int_mb(i_indx+2*ii), 859 > dbl_mb(i_coef+ii), 860 > K0,gcut0,gamma0,dbl_mb(i_c),gamm,enb,gx) 861 if(oprint) then 862 write(luout,*) 863 write(luout,'(A,I4)') "Adding bondings penalty # ",i 864 write(luout,'(A,F18.9,F18.9,E14.6)') 865 > " - spring parameters (K0,gcut0,gamma0):", 866 > K0,gcut0,gamma0 867 write(luout,'(A,E18.6)') 868 > " - gamma :",gamm 869 write(luout,'(A,F18.9)') 870 > " - penalty energy :",enb 871 write(luout,'(A,F18.9)') 872 > " - previous/base energy :",energy 873 write(luout,'(A,F18.9)') 874 > " - penalized energy :",energy+enb 875 write(luout,*) " - coefficient index1 index2" 876 do j=1,n 877 write(luout,'(F16.6,I7,I7)') dbl_mb(i_coef+ii+j-1), 878 > int_mb(i_indx+2*ii+2*(j-1)), 879 > int_mb(i_indx+2*ii+2*(j-1)+1) 880 end do 881 write(luout,*) 882 end if 883 884 ii = ii + n 885 e=e+enb 886 end do 887 energy = energy + e 888c 889 call util_print_pop() 890 return 891 end 892 893 894 895 896 897 898c 899 subroutine cons_add_hdihed_egrad(rtdb,energy,gx) 900 implicit none 901#include "errquit.fh" 902#include "global.fh" 903#include "mafdecls.fh" 904#include "rtdb.fh" 905#include "util.fh" 906#include "stdio.fh" 907#include "cons.fh" 908#include "cons_data.fh" 909c 910 integer rtdb 911 double precision energy 912 double precision gx(*) 913c 914 logical status 915 logical oprint 916 integer inb 917 integer iat,jat,kat,lat 918 double precision r 919 double precision r0,k 920 double precision e,enb 921 double precision f1(3),f2(3),f3(3),f4(3) 922 integer i_c,i_r,i_k,i_id 923 integer nc 924 integer idm 925 integer i 926c 927 idm = 4 928c 929 call util_print_push() 930 call util_print_rtdb_load(rtdb,'cons') 931 oprint = util_print('energy', print_default) 932 oprint =oprint .and. (ga_nodeid().eq.0) 933c 934 call cons_get_hdihed_n(nc) 935 i_c = cons_get_i_c() 936 i_id = cons_get_i_hdihed_id() 937 i_k = cons_get_i_hdihed_k() 938 i_r = cons_get_i_hdihed_r() 939 940 if(oprint) then 941 write(luout,*)"Processing dihedral constraints" 942 write(luout,1004) 943 end if 944 do inb=1,nc 945 iat=int_mb(i_id+idm*(inb-1)) 946 jat=int_mb(i_id+idm*(inb-1)+1) 947 kat=int_mb(i_id+idm*(inb-1)+2) 948 lat=int_mb(i_id+idm*(inb-1)+3) 949 r0 =dbl_mb(i_r+inb-1) 950 k =dbl_mb(i_k+inb-1) 951 952 call cons_dihed_force(k,r0, 953 > dbl_mb(i_c+(iat-1)*3), 954 > dbl_mb(i_c+(jat-1)*3), 955 > dbl_mb(i_c+(kat-1)*3), 956 > dbl_mb(i_c+(lat-1)*3), 957 > r, 958 > enb, 959 > f1, 960 > f2, 961 > f3, 962 > f4) 963 964 if(oprint) then 965 write(luout,1005) iat,jat,kat,lat,k, 966 > r0*rad_to_deg,r*rad_to_deg, 967 > enb,f1(1),f2(1),f3(1),f4(1) 968 write(luout,1006) f1(2),f2(2),f3(2),f4(2) 969 write(luout,1006) f1(3),f2(3),f3(3),f4(3) 970c write(6,1002) "dihedral forces (i,j,k,l):" 971c write(6,1001) "f_i", iat,(f1(i),i=1,3) 972c write(6,1001) "f_j", jat,(f2(i),i=1,3) 973c write(6,1001) "f_k", kat,(f3(i),i=1,3) 974c write(6,1001) "f_l", lat,(f4(i),i=1,3) 975 end if 976 do i=1,3 977 gx((iat-1)*3+i)=gx((iat-1)*3+i)+f1(i) 978 gx((jat-1)*3+i)=gx((jat-1)*3+i)+f2(i) 979 gx((kat-1)*3+i)=gx((kat-1)*3+i)+f3(i) 980 gx((lat-1)*3+i)=gx((lat-1)*3+i)+f4(i) 981 end do 982 energy = energy + enb 983 end do 984c 985 call util_print_pop() 9861002 FORMAT(1X,A,/,T6,"atom",T13,"fx",T24,"fy",T35,"fz") 9871001 FORMAT(1X,A,1X,I4,3(1X,F10.6)) 9881004 FORMAT(T5,"i",4X,"j",4X,"k",4X,"l",T23,"Kphi",T31,"phi0", 989 > T39,"phi",T47,"Energy",T55,"f1",T62,"f2",T69"f3",T76,"f4",/, 990 > T5,76("_")) 9911005 FORMAT( 4(1X,I4),3(1X,F7.3),1X,5(F7.3)) 9921006 FORMAT( T53,4F7.3) 993 return 994 end 995c 996 subroutine cons_add_hdihed_energy(rtdb,energy) 997 implicit none 998#include "errquit.fh" 999#include "global.fh" 1000#include "mafdecls.fh" 1001#include "rtdb.fh" 1002#include "util.fh" 1003#include "stdio.fh" 1004#include "cons.fh" 1005#include "cons_data.fh" 1006c 1007 integer rtdb 1008 double precision energy 1009 logical status 1010 logical oprint 1011 integer inb 1012 integer iat,jat,kat,lat 1013 double precision r 1014 double precision r0,k 1015 double precision e,enb 1016 integer i_c,i_r,i_k,i_id 1017 integer nc 1018 integer idm 1019c 1020 idm = 4 1021c 1022 call util_print_push() 1023 call util_print_rtdb_load(rtdb,'cons') 1024 oprint = util_print('energy', print_default) 1025 oprint =oprint .and. (ga_nodeid().eq.0) 1026c 1027 call cons_get_hdihed_n(nc) 1028 i_c = cons_get_i_c() 1029 i_id = cons_get_i_hdihed_id() 1030 i_k = cons_get_i_hdihed_k() 1031 i_r = cons_get_i_hdihed_r() 1032 1033 if(oprint) then 1034 write(luout,*)"Processing dihedral constraints" 1035 write(luout,1004) 1036 end if 1037 do inb=1,nc 1038 iat=int_mb(i_id+idm*(inb-1)) 1039 jat=int_mb(i_id+idm*(inb-1)+1) 1040 kat=int_mb(i_id+idm*(inb-1)+2) 1041 lat=int_mb(i_id+idm*(inb-1)+3) 1042 r0 =dbl_mb(i_r+inb-1) 1043 k =dbl_mb(i_k+inb-1) 1044 1045 call cons_dihed_energy(k,r0, 1046 > dbl_mb(i_c+(iat-1)*3), 1047 > dbl_mb(i_c+(jat-1)*3), 1048 > dbl_mb(i_c+(kat-1)*3), 1049 > dbl_mb(i_c+(lat-1)*3), 1050 > r, 1051 > enb) 1052 1053 if(oprint) then 1054 write(luout,1005) iat,jat,kat,lat,k, 1055 > r0*rad_to_deg,r*rad_to_deg, 1056 > enb 1057 end if 1058 1059 energy=energy+enb 1060 end do 1061c 1062 call util_print_pop() 1063 return 10641004 FORMAT(T5,"i",4X,"j",4X,"k",4X,"l",T23,"Kphi",T31,"phi0", 1065 > T39,"phi",T47,"Energy",/, 1066 > T5,49("_")) 10671005 FORMAT( 4(1X,I4),3(1X,F7.3),1X,G12.6) 1068 end 1069c 1070 subroutine cons_add_hbond_egrad(rtdb,energy,gx) 1071 implicit none 1072#include "errquit.fh" 1073#include "mafdecls.fh" 1074#include "global.fh" 1075#include "rtdb.fh" 1076#include "stdio.fh" 1077#include "util.fh" 1078#include "cons.fh" 1079c 1080 integer rtdb 1081 double precision gx(*) 1082 double precision energy 1083c local variables 1084c --------------- 1085 logical status 1086 logical oprint 1087 integer i 1088 integer inb,nhb 1089 integer iat,jat 1090 double precision r 1091 double precision r0,k 1092 double precision e,f(3) 1093 integer i_c,i_rhb,i_khb,i_ijhb 1094c 1095 call util_print_push() 1096 call util_print_rtdb_load(rtdb,'cons') 1097 oprint = util_print('energy', print_default) 1098 oprint =oprint .and. (ga_nodeid().eq.0) 1099c 1100 call cons_get_hbond_nhb(nhb) 1101 i_c = cons_get_i_c() 1102 i_ijhb = cons_get_i_hbond_id() 1103 i_khb = cons_get_i_hbond_k() 1104 i_rhb = cons_get_i_hbond_r() 1105 1106c 1107 e=0.0d0 1108 do inb=1,nhb 1109 iat=int_mb(i_ijhb+2*(inb-1)) 1110 jat=int_mb(i_ijhb+2*(inb-1)+1) 1111 r0 =dbl_mb(i_rhb+inb-1) 1112 k =dbl_mb(i_khb+inb-1) 1113c 1114 call cons_spring_force(k,r0, 1115 > dbl_mb(i_c+(iat-1)*3), 1116 > dbl_mb(i_c+(jat-1)*3), 1117 > r,e,f) 1118 1119 if(oprint) then 1120 write(luout,*)"adding spring # ",inb 1121 write(luout,*)" spring parameters (i,j,k,r0):",iat,jat,k,r0 1122 write(luout,*)" spring length and energy :",r,e 1123 write(luout,*)" spring forces :",(f(i),i=1,3) 1124 write(luout,*)" spring deriv :", 1125 > -2.0*k*(r-r0) 1126 end if 1127 1128 if(rtdb_get(rtdb, 'cons:simulate', mt_log, 1, status)) then 1129 do i=1,3 1130 gx((iat-1)*3+i)=f(i) 1131 gx((jat-1)*3+i)=-f(i) 1132 end do 1133 energy = e 1134 else 1135 do i=1,3 1136 gx((iat-1)*3+i)=gx((iat-1)*3+i)+f(i) 1137 gx((jat-1)*3+i)=gx((jat-1)*3+i)-f(i) 1138 end do 1139 energy = energy+e 1140 end if 1141 1142 end do 1143 call util_print_pop() 1144 1145 return 1146 end 1147c 1148 subroutine cons_spring_force(k,r0,r1,r2,r,energy,f) 1149 implicit none 1150#include "errquit.fh" 1151#include "mafdecls.fh" 1152#include "rtdb.fh" 1153#include "util.fh" 1154#include "cons_data.fh" 1155c 1156 double precision k,r0 1157 double precision r1(3) 1158 double precision r2(3) 1159 double precision r 1160 double precision energy 1161 double precision f(3) 1162c 1163 integer i 1164c 1165 r=(r1(1)-r2(1))**2+ 1166 > (r1(2)-r2(2))**2+ 1167 > (r1(3)-r2(3))**2 1168 r=sqrt(r) 1169 energy=k*(r-r0)**2 1170 1171 do i=1,3 1172 f(i)=2*k*(r-r0)* 1173 & (r1(i)-r2(i))/r 1174 end do 1175 1176 return 1177 end 1178 1179 subroutine cons_spring_energy(k,r0,r1,r2,r,energy) 1180 implicit none 1181#include "errquit.fh" 1182#include "mafdecls.fh" 1183#include "rtdb.fh" 1184#include "util.fh" 1185#include "cons_data.fh" 1186c 1187 double precision k,r0 1188 double precision r1(3) 1189 double precision r2(3) 1190 double precision r 1191 double precision energy 1192c 1193 integer i 1194c 1195 r=(r1(1)-r2(1))**2+ 1196 > (r1(2)-r2(2))**2+ 1197 > (r1(3)-r2(3))**2 1198 r=sqrt(r) 1199 energy=k*(r-r0)**2 1200 return 1201 end 1202 1203 1204 subroutine cons_bondings_energy(n,indx,coef,K0,gamma0,rion, 1205 > gamm,energy) 1206 implicit none 1207#include "errquit.fh" 1208#include "mafdecls.fh" 1209#include "rtdb.fh" 1210#include "util.fh" 1211#include "cons_data.fh" 1212c 1213 integer n,indx(*) 1214 double precision coef(*),K0,gamma0 1215 double precision rion(3,*) 1216 double precision gamm,energy 1217c 1218 integer i,i0,i1 1219 double precision x,y,z,dist 1220c 1221 gamm = 0.0d0 1222 do i=1,n 1223 i0 = indx(2*i-1) 1224 i1 = indx(2*i) 1225 x = rion(1,i0) - rion(1,i1) 1226 y = rion(2,i0) - rion(2,i1) 1227 z = rion(3,i0) - rion(3,i1) 1228 dist = dsqrt(x*x + y*y + z*z) 1229 gamm = gamm + coef(i)*dist 1230 end do 1231 energy = K0*(gamm-gamma0)**2 1232 1233 return 1234 end 1235 1236 subroutine cons_bondings_force(n,indx,coef,K0,gamma0,rion, 1237 > gamm,energy,f) 1238 implicit none 1239#include "errquit.fh" 1240#include "mafdecls.fh" 1241#include "rtdb.fh" 1242#include "util.fh" 1243#include "cons_data.fh" 1244c 1245 integer n,indx(*) 1246 double precision coef(*),K0,gamma0 1247 double precision rion(3,*) 1248 double precision gamm,energy,f(3,*) 1249c 1250 integer i,ii,i0,i1 1251 double precision x,y,z,dist,denergy 1252c 1253 gamm = 0.0d0 1254 do i=1,n 1255 i0 = indx(2*i-1) 1256 i1 = indx(2*i) 1257 x = rion(1,i0) - rion(1,i1) 1258 y = rion(2,i0) - rion(2,i1) 1259 z = rion(3,i0) - rion(3,i1) 1260 dist = dsqrt(x*x + y*y + z*z) 1261 gamm = gamm + coef(i)*dist 1262 end do 1263 energy = K0*(gamm-gamma0)**2 1264 denergy = 2.0d0*K0*(gamm-gamma0) 1265 1266 do i=1,n 1267 i0 = indx(2*i-1) 1268 i1 = indx(2*i) 1269 x = rion(1,i0) - rion(1,i1) 1270 y = rion(2,i0) - rion(2,i1) 1271 z = rion(3,i0) - rion(3,i1) 1272 dist = dsqrt(x*x + y*y + z*z) 1273 f(1,i0) = f(1,i0) + denergy*coef(i)*(x/dist) 1274 f(2,i0) = f(2,i0) + denergy*coef(i)*(y/dist) 1275 f(3,i0) = f(3,i0) + denergy*coef(i)*(z/dist) 1276 f(1,i1) = f(1,i1) - denergy*coef(i)*(x/dist) 1277 f(2,i1) = f(2,i1) - denergy*coef(i)*(y/dist) 1278 f(3,i1) = f(3,i1) - denergy*coef(i)*(z/dist) 1279 end do 1280 1281 return 1282 end 1283 1284 1285 1286 subroutine cons_pbondings_energy(n,indx,coef,K0,gcut0,gamma0,rion, 1287 > gamm,energy) 1288 implicit none 1289#include "errquit.fh" 1290#include "mafdecls.fh" 1291#include "rtdb.fh" 1292#include "util.fh" 1293#include "cons_data.fh" 1294c 1295 integer n,indx(*) 1296 double precision coef(*),K0,gcut0,gamma0 1297 double precision rion(3,*) 1298 double precision gamm,energy 1299c 1300 integer i,i0,i1 1301 double precision x,y,z,dist 1302c 1303 gamm = 0.0d0 1304 do i=1,n 1305 i0 = indx(2*i-1) 1306 i1 = indx(2*i) 1307 x = rion(1,i0) - rion(1,i1) 1308 y = rion(2,i0) - rion(2,i1) 1309 z = rion(3,i0) - rion(3,i1) 1310 dist = dsqrt(x*x + y*y + z*z) 1311 gamm = gamm + coef(i)*dist 1312 end do 1313 energy = K0*exp(-((gamm-gamma0)/gcut0)**2) 1314 1315 return 1316 end 1317 1318 1319 1320 subroutine cons_pbondings_force(n,indx,coef,K0,gcut0,gamma0,rion, 1321 > gamm,energy,f) 1322 implicit none 1323#include "errquit.fh" 1324#include "mafdecls.fh" 1325#include "rtdb.fh" 1326#include "util.fh" 1327#include "cons_data.fh" 1328c 1329 integer n,indx(*) 1330 double precision coef(*),K0,gcut0,gamma0 1331 double precision rion(3,*) 1332 double precision gamm,energy,f(3,*) 1333c 1334 integer i,ii,i0,i1 1335 double precision x,y,z,dist,denergy 1336c 1337 gamm = 0.0d0 1338 do i=1,n 1339 i0 = indx(2*i-1) 1340 i1 = indx(2*i) 1341 x = rion(1,i0) - rion(1,i1) 1342 y = rion(2,i0) - rion(2,i1) 1343 z = rion(3,i0) - rion(3,i1) 1344 dist = dsqrt(x*x + y*y + z*z) 1345 gamm = gamm + coef(i)*dist 1346 end do 1347 energy = K0*exp(-((gamm-gamma0)/gcut0)**2) 1348 !denergy = 2.0d0*K0*(gamm-gamma0) 1349 denergy = -2.0d0*K0*(gamm-gamma0)*exp(-((gamm-gamma0)/gcut0)**2) 1350 > / gcut0**2 1351 1352 do i=1,n 1353 i0 = indx(2*i-1) 1354 i1 = indx(2*i) 1355 x = rion(1,i0) - rion(1,i1) 1356 y = rion(2,i0) - rion(2,i1) 1357 z = rion(3,i0) - rion(3,i1) 1358 dist = dsqrt(x*x + y*y + z*z) 1359 f(1,i0) = f(1,i0) + denergy*coef(i)*(x/dist) 1360 f(2,i0) = f(2,i0) + denergy*coef(i)*(y/dist) 1361 f(3,i0) = f(3,i0) + denergy*coef(i)*(z/dist) 1362 f(1,i1) = f(1,i1) - denergy*coef(i)*(x/dist) 1363 f(2,i1) = f(2,i1) - denergy*coef(i)*(y/dist) 1364 f(3,i1) = f(3,i1) - denergy*coef(i)*(z/dist) 1365 end do 1366 1367 return 1368 end 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 subroutine cons_dihed_energy(k,phi0,r1,r2,r3,r4,phi,energy) 1380 implicit none 1381#include "errquit.fh" 1382#include "mafdecls.fh" 1383#include "rtdb.fh" 1384#include "util.fh" 1385#include "cons_data.fh" 1386c 1387 double precision k,phi0 1388 double precision r1(3) 1389 double precision r2(3) 1390 double precision r3(3) 1391 double precision r4(3) 1392 double precision phi 1393 double precision energy 1394c 1395 integer i 1396c 1397c calculate dihedral angle 1398c ------------------------ 1399 call cons_dihed(r1,r2,r3,r4,phi,"rads") 1400 energy = 0.5d0*k*(phi-phi0)**2 1401 return 1402 end 1403 1404 subroutine cons_dihed_force(k,phi0,r1,r2,r3,r4, 1405 > phi,energy,f1,f2,f3,f4) 1406 implicit none 1407#include "errquit.fh" 1408#include "mafdecls.fh" 1409#include "rtdb.fh" 1410#include "util.fh" 1411#include "cons_data.fh" 1412c 1413 double precision k,phi0 1414 double precision r1(3) 1415 double precision r2(3) 1416 double precision r3(3) 1417 double precision r4(3) 1418 double precision phi 1419 double precision energy 1420 double precision f1(3) 1421 double precision f2(3) 1422 double precision f3(3) 1423 double precision f4(3) 1424c 1425 integer i 1426 double precision a 1427c 1428c calculate dihedral angle 1429c ------------------------ 1430 call cons_dihed(r1,r2,r3,r4,phi,"rads") 1431 energy = 0.5d0*k*(phi-phi0)**2 1432 call cons_dihed_deriv(r1,r2,r3,r4,f1,f2,f3,f4,"rads") 1433 a = k*(phi-phi0) 1434 do i=1,3 1435 f1(i) = f1(i)*a 1436 f2(i) = f2(i)*a 1437 f3(i) = f3(i)*a 1438 f4(i) = f4(i)*a 1439 end do 1440 return 1441 end 1442 1443 subroutine cons_load_hdihed(namespace,rtdb) 1444 implicit none 1445#include "errquit.fh" 1446#include "mafdecls.fh" 1447#include "rtdb.fh" 1448#include "util.fh" 1449#include "geom.fh" 1450c 1451 character*(*) namespace 1452 integer rtdb 1453c 1454 logical status 1455 integer nhd 1456 integer i_rhd,i_khd,i_ijhd 1457 integer h_rhd,h_khd,h_ijhd 1458 character*255 tag_id 1459 character*255 tag_n 1460 character*255 tag_r 1461 character*255 tag_k 1462 integer ma_type,ma_n 1463c 1464 call cons_unload_hdihed() 1465c 1466 call cons_hdihed_id_tag(namespace,tag_id) 1467 call cons_hdihed_n_tag(namespace,tag_n) 1468 call cons_hdihed_k_tag(namespace,tag_k) 1469 call cons_hdihed_r_tag(namespace,tag_r) 1470c 1471c load harmonic constraints 1472c 1473 status=rtdb_get(rtdb,tag_n,mt_int,1,nhd) 1474c 1475c exit if no harm dihedrals 1476 if(.not.status .or. nhd.eq.0) then 1477c call errquit( 1478c > 'cons_load_hdiheds: unable to get number of harm bonds:', 1479c > nhd, MA_ERR) 1480 return 1481 end if 1482 1483 if ( .not. rtdb_ma_get( rtdb, tag_id,ma_type, ma_n, 1484 & h_ijhd) ) call errquit( 1485 & 'cons_load_hdiheds: unable to allocate cons scratch space', 1486 & 2*nhd, MA_ERR) 1487 if ( .not. ma_get_index(h_ijhd, 1488 & i_ijhd) ) call errquit( 1489 & 'cons_load_hdiheds: unable to allocate cons scratch space', 1490 & 2*nhd, MA_ERR) 1491 1492 if ( .not. rtdb_ma_get( rtdb, tag_k,ma_type, ma_n, 1493 & h_khd) ) call errquit( 1494 & 'cons_load_hdiheds: unable to allocate cons scratch space', 1495 & 2*nhd, MA_ERR) 1496 if ( .not. ma_get_index(h_khd, 1497 & i_khd) ) call errquit( 1498 & 'cons_load_hdiheds: unable to allocate cons scratch space', 1499 & 2*nhd, MA_ERR) 1500 1501 if ( .not. rtdb_ma_get( rtdb, tag_r,ma_type, ma_n, 1502 & h_rhd) ) call errquit( 1503 & 'cons_load_hdiheds: unable to allocate cons scratch space', 1504 & 2*nhd, MA_ERR) 1505 if ( .not. ma_get_index(h_rhd, 1506 & i_rhd) ) call errquit( 1507 & 'cons_load_hdiheds: unable to allocate cons scratch space', 1508 & 2*nhd, MA_ERR) 1509 1510 1511 call cons_set_hdihed_n(nhd) 1512 call cons_set_h_hdihed_id(h_ijhd) 1513 call cons_set_h_hdihed_k(h_khd) 1514 call cons_set_h_hdihed_r(h_rhd) 1515 1516 end 1517c 1518 subroutine cons_unload_hdihed() 1519 implicit none 1520#include "errquit.fh" 1521#include "mafdecls.fh" 1522#include "cons.fh" 1523 integer nhd,h_rhd,h_khd,h_idhd 1524 1525 h_idhd = cons_get_h_hdihed_id() 1526 h_khd = cons_get_h_hdihed_k() 1527 h_rhd = cons_get_h_hdihed_r() 1528c 1529c unload harmonic constraints if any 1530c 1531 call cons_get_hdihed_n(nhd) 1532 if(nhd.gt.0) then 1533 if (.not. ma_free_heap(h_rhd) ) call errquit( 1534 & 'cons_bond_input: unable to free h_rhd', 1535 & 0, MA_ERR) 1536 if (.not. ma_free_heap(h_khd) ) call errquit( 1537 & 'cons_bond_input: unable to free h_khd', 1538 & 0, MA_ERR) 1539 if (.not. ma_free_heap(h_idhd) ) call errquit( 1540 & 'cons_bond_input: unable to free h_ijhd', 1541 & 0, MA_ERR) 1542 call cons_set_hdihed_n(0) 1543 end if 1544 end 1545c $Id$ 1546