1********************************************************************** 2* * 3* REAXFF Reactive force field program * 4* * 5* Developed and written by Adri van Duin, duin@wag.caltech.edu * 6* * 7* Copyright (c) 2001-2010 California Institute of Technology * 8* * 9* This is an open-source program. Feel free to modify its * 10* contents. Please keep me informed of any useful modification * 11* or addition that you made. Please do not distribute this * 12* program to others; if people are interested in obtaining * 13* a copy of this program let them contact me first. * 14* * 15********************************************************************** 16******************************************************************** 17 18 subroutine calval 19 20********************************************************************** 21#include "cbka.blk" 22#include "cbkc.blk" 23#include "cbkdhdc.blk" 24#include "cbkdrdc.blk" 25#include "cbkh.blk" 26#include "cbkrbo.blk" 27#include "cbkvalence.blk" 28#include "cellcoord.blk" 29#include "control.blk" 30 dimension a(3),b(3),j(3),dradc(3,3),drbdc(3,3),dtdc(3,3), 31 $dargdc(3,3),dndc(3,3),dadc(3),dbdc(3) 32********************************************************************** 33* * 34* Calculate valency angles and their derivatives to cartesian * 35* coordinates * 36* Valency angle energies are calculated in valang * 37* * 38********************************************************************** 39********************************************************************** 40* Description of variables used in this routine. 41* 42* ndebug: stored in cbka.blk; control-parameter 43* third: local variable 44* twothird: local variable 45* dadc(3): local array; stores derivative distance to cartesians 46* dbdc(3): local array; stores derivative distance to cartesians 47* i1: local do-loop counter 48* i2: local do-loop counter 49* k1: local do-loop counter 50* k2: local do-loop counter 51* dradc(3,3): local array; stores derivatives bond lengths to 52* cartesians 53* drbdc(3,3): local array; stores derivatives bond lengths to 54* cartesians 55* nval: stored in cbka.blk; number of valence angles 56* ity: local integer; atom type 57* iv(nvalmax,6): stored in cbka.blk; valence angle identifiers 58* j(3): local integer array; stores valence angle atom numbers 59* la: local integer: stores bond numbers in valence angle 60* lb: local integer: stores bond numbers in valence angle 61* ivl1: local integer; stores symmetric copy number of bond 62* ivl2: local integer; stores symmetric copy number of bond 63* ibsym(nbomax): stored in cbka.blk; symmetric copy number of bond 64* isign1: local integer; -1 or 1 65* isign2: local integer; -1 or 1 66* rla: local variable; stores bond length for bond la 67* rlb: local variable; stores bond length for bond lb 68* rbo(nbomax): stored in cbka.blk; stores bond lengths 69* ix1,iy1,iz1,ix2,iy2,iz2: local integers; periodic cell shifts 70* a(3): local variable; distance in x,y and z-direction between atoms 71* b(3): local variable; distance in x,y and z-direction between atoms 72* c(nat,3): stored in cbka.blk; cartesian coordinate array 73* tm11,tm21,tm22,tm31,tm32,tm33: stored in cbka.blk; periodic cell 74* matrix 75* poem: local variable; product of bond lengths 76* tel: local variable; cross-product of x,y and z-interatomic 77* distances 78* arg: local variable; cosine of angle between bonds a and b 79* arg2: local variable; square of arg 80* s1ma22: local variable; used to check whether angle gets to 180 81* degrees 82* s1ma2: local variable; square root of s1ma22 83* hl: local variable; angle (in radians) between bonds a and b 84* h(nvamax): stored in cbka.blk; angle (in radians) between bonds a 85* and b 86* ib(nbomax,3): stored in cbka.blk: bond distance identifiers 87* drdc(nbomax,3,2): stored in cbka.blk; derivatives bond distances 88* to cartesian coordinates 89* dndc(3,3): local variable; temporary storage for calculating 90* derivatives of valence angle to cartesians 91* dtdc(3,3): local variable; temporary storage for calculating 92* derivatives of valence angle to cartesians 93* dargdc(3,3): local variable; temporary storage for calculating 94* derivatives of valence angle to cartesians 95* dhdc(nvamax,3,3): stored in cbka.blk; derivatives of valence angle 96* to cartesians 97* 98********************************************************************** 99c$$$ if (ndebug.eq.1) then 100c$$$C open (65,file='fort.65',status='unknown',access='append') 101c$$$ write (65,*) 'In calval' 102c$$$ call timer(65) 103c$$$ close (65) 104c$$$ end if 105 106 third=1.0/3.0 107 twothird=2.0/3.0 108 dadc(1)=-1.0 109 dadc(2)=1.0 110 dadc(3)=0.0 111 dbdc(1)=0.0 112 dbdc(2)=1.0 113 dbdc(3)=-1.0 114 do k1=1,3 115 do k2=1,3 116 dradc(k1,k2)=0.0 117 drbdc(k1,k2)=0.0 118 end do 119 end do 120 if (nval.eq.0) return 121 122 do 10 i1=1,nval 123 ity=iv(i1,1) 124 j(1)=iv(i1,2) 125 j(2)=iv(i1,3) 126 j(3)=iv(i1,4) 127********************************************************************** 128* * 129* Determine valency angle * 130* * 131********************************************************************** 132 la=iv(i1,5) 133 lb=iv(i1,6) 134 ivl1=ibsym(la) 135 ivl2=ibsym(lb) 136 isign1=1 137 isign2=1 138 rla=rbo(la) 139 rlb=rbo(lb) 140 141 call dista2(j(2),j(1),dis,a(1),a(2),a(3)) 142 call dista2(j(2),j(3),dis,b(1),b(2),b(3)) 143 144 poem=rla*rlb 145 tel=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) 146 arg=tel/poem 147 arg2=arg*arg 148 s1ma22=1.0-arg2 149 if (s1ma22.lt.1.0d-10) s1ma22=1.0d-10 150 s1ma2=sqrt(s1ma22) 151 if (arg.gt.1.0) arg=1.0 152 if (arg.lt.-1.0) arg=-1.0 153 hl=acos(arg) 154 h(i1)=hl 155********************************************************************** 156* * 157* Calculate derivative valency angle to cartesian coordinates * 158* * 159********************************************************************** 160 if (j(1).eq.ib(la,2)) then 161 do k1=1,3 162 dradc(k1,1)=drdc(la,k1,1) 163 dradc(k1,2)=drdc(la,k1,2) 164 end do 165 else 166 do k1=1,3 167 dradc(k1,1)=drdc(la,k1,2) 168 dradc(k1,2)=drdc(la,k1,1) 169 end do 170 end if 171 if (j(2).eq.ib(lb,2)) then 172 do k1=1,3 173 drbdc(k1,2)=drdc(lb,k1,1) 174 drbdc(k1,3)=drdc(lb,k1,2) 175 end do 176 else 177 do k1=1,3 178 drbdc(k1,2)=drdc(lb,k1,2) 179 drbdc(k1,3)=drdc(lb,k1,1) 180 end do 181 end if 182 do k1=1,3 183 do k2=1,3 184 dndc(k1,k2)=rla*drbdc(k1,k2)+rlb*dradc(k1,k2) 185 dtdc(k1,k2)=a(k1)*dbdc(k2)+b(k1)*dadc(k2) 186 dargdc(k1,k2)=(dtdc(k1,k2)-arg*dndc(k1,k2))/poem 187 dhdc(i1,k1,k2)=-dargdc(k1,k2)/s1ma2 188 end do 189 end do 190 191 10 continue 192 193 return 194 end 195********************************************************************** 196********************************************************************** 197 198 subroutine boncor 199 200********************************************************************** 201#include "cbka.blk" 202#include "cbkabo.blk" 203#include "cbkc.blk" 204#include "cbkbo.blk" 205#include "cbkboncor.blk" 206#include "cbkbosi.blk" 207#include "cbkbopi.blk" 208#include "cbkbopi2.blk" 209#include "cbkconst.blk" 210#include "cbkdbopi2ndc.blk" 211#include "cbkdbopidc.blk" 212#include "cbkdbopindc.blk" 213#include "cbkff.blk" 214#include "cbkia.blk" 215#include "cbkidbo.blk" 216#include "cbknubon2.blk" 217#include "cbkrbo.blk" 218#include "control.blk" 219#include "small.blk" 220#include "cbkdbodc.blk" 221c$$$ if (ndebug.eq.1) then 222c$$$C open (65,file='fort.65',status='unknown',access='append') 223c$$$ write (65,*) 'In boncor' 224c$$$ call timer(65) 225c$$$ close (65) 226c$$$ end if 227********************************************************************** 228* * 229* Correction for overcoordination and 1-3 bond orders * 230* * 231********************************************************************** 232********************************************************************** 233* Description of variables used in this routine. 234* 235* ndebug: stored in cbka.blk; control-parameter 236* i1: local do-loop counter 237* i2: local do-loop counter 238* k1: local do-loop counter 239* k2: local do-loop counter 240* nbon: stored in cbka.blk; number of bonds in system 241* ibt: local integer; stores bond type 242* ib(nbomax,3): stored in cbka.blk: bond distance identifiers 243* j1: local integer; stores atom number 1st atom in bond 244* j2: local integer; stores atom number 2nd atom in bond 245* ovc(nbotym): stored in cbka.blk: force field parameter for 246* overcoordination correction 247* v13cor(nbotym): stored in cbka.blk: force field parameter for 248* 1-3 bond order correction 249* idbo1(nbomax): stored in cbka.blk; number of atoms in the 250* derivative of the bond order 251* idbo(nbomax,2*mbond): stored in cbka.blk; atom numbers of the 252* atoms in the derivative of the bond order 253* dbondc(nbomax,3,2*mbond): stored in cbka.blk; derivative of 254* corrected total bond orders to cartesians 255* dbosindc(nbomax,3,2*mbond): stored in cbka.blk; derivative of 256* corrected sigma bond orders to cartesians 257* dbopindc(nbomax,3,2*mbond): stored in cbka.blk; derivative of 258* corrected pi bond orders to cartesians 259* dbopi2ndc(nbomax,3,2*mbond): stored in cbka.blk; derivative of 260* corrected double pi bond orders to cartesians 261* dbodc(nbomax,3,2): stored in cbka.blk; derivative of 262* uncorrected total bond orders to cartesians 263* dbosidc(nbomax,3,2): stored in cbka.blk; derivative of 264* uncorrected sigma bond orders to cartesians 265* dbopidc(nbomax,3,2): stored in cbka.blk; derivative of 266* uncorrected pi bond orders to cartesians 267* dbopi2dc(nbomax,3,2): stored in cbka.blk; derivative of 268* uncorrected double pi bond orders to cartesians 269* boo: local variable; storage of uncorrected total bond order 270* bo(nbomax): stored in cbka.blk; total bond order 271* bopi(nbomax): stored in cbka.blk; pi bond order 272* bopi2(nbomax): stored in cbka.blk; double pi bond order 273* bopio: local variable; storage of uncorrected pi bond order 274* bopi2o: local variable; storage of uncorrected double pi bond order 275* iti: local integer; atom type first atom in bond 276* itj: local integer; atom type second atom in bond 277* ia(nat,mbond+3): stored in cbka.blk; connection table without bond 278* order cutoff 279* aboi: local variable: total bond order around atom i 280* aboj: local variable: total bond order around atom j 281* abo(nat): stored in cbka.blk; total bond order around atoms 282* vp131: local variable; force field cross-term 283* vp132: local variable; force field cross-term 284* vp133: local variable; force field cross-term 285* bo131(nsort): stored in cbka.blk; force field parameter for 1-3 286* bond order correction 287* bo132(nsort): stored in cbka.blk; force field parameter for 1-3 288* bond order correction 289* bo133(nsort): stored in cbka.blk; force field parameter for 1-3 290* bond order correction 291* corrtot:local variable; total correction on bond order 292* dbodsboi1: local variable; derivative of bond order to sum of bond 293* orders around atom i 294* dbodsboj1: local variable; derivative of bond order to sum of bond 295* orders around atom j 296* ovi: local variable; overcoordination on atom i 297* ovj: local variable; overcoordination on atom j 298* aval(nat): stored in cbka.blk; nr. of valence electrons on atom 299* exphu1: local variable; stores exponential 300* exphu2: local variable; stores exponential 301* exp11: local variable; stores exponential 302* exp21: local variable; stores exponential 303* vpar(npamax): stored in cbka.blk: general parameters 304* exphu12: local variable; stores sum of exponential 305* ovcor: local variable; temporary storage for BO/ovcor corr. 306* huli: local variable; temporary storage for BO/ovcor corr. 307* hulj: local variable; temporary storage for BO/ovcor corr. 308* corr1: local variable; temporary storage for BO/ovcor corr. 309* corr2: local variable; temporary storage for BO/ovcor corr. 310* dbodsboi2: local variable; derivative of 1-3 BO correction to sum 311* of bond orders around atom i 312* dbodsboj2: local variable; derivative of 1-3 BO correction to sum 313* of bond orders around atom i 314* bocor1: local variable; 1-3 bond order correction 315* bocor2: local variable; 1-3 bond order correction 316* ovi2: local variable; overcoordination on atom i with reference to 317* total number of electrons on atom i, including lone 318* pairs 319* ovj2: local variable; overcoordination on atom j with reference to 320* total number of electrons on atom j, including lone 321* pairs 322* valf(nsort): stored in cbka.blk; total number of electrons on 323* atom, including lone pairs 324* cor1: local variable; temporary storage for BO/1-3 bond corr. 325* cor2: local variable; temporary storage for BO/1-3 bond corr. 326* exphu3: local variable; storage exponential 327* exphu4: local variable; storage exponential 328* corrtot2: local variable; square of corrtot 329* dbodboo: local variable; derivative of corrected total bond order to 330* uncorrected bond order 331* dbopidbopio: local variable; derivative of corrected pi bond order 332* to uncorrected pi bond order 333* dbopidboo: local variable; derivative of corrected pi bond order 334* to uncorrected total bond order 335* dbopi2dbopi2o: local variable; derivative of corrected double pi bond order 336* to uncorrected double pi bond order 337* dbopi2dboo: local variable; derivative of corrected double pi bond order 338* to uncorrected total bond order 339* dbodsboit: local variable; derivative of total bond order to sum 340* of bond orders around atom i 341* dbodsbojt: local variable; derivative of total bond order to sum 342* of bond orders around atom j 343* vhui: local variable; temporary storage 344* vhuj: local variable; temporary storage 345* dbopidsboit: local variable; derivative of pi bond order to sum 346* of bond orders around atom i 347* dbopidsbojt: local variable; derivative of pi bond order to sum 348* of bond orders around atom j 349* dbopi2dsboit: local variable; derivative of pi bond order to sum 350* of bond orders around atom i 351* dbopi2dsbojt: local variable; derivative of pi bond order to sum 352* of bond orders around atom j 353* nco: local integer; counter for number of atoms in derivative 354* ihl: local integer; helps to access right dbodc-term 355* nubon2(nat,mbond): stored in cbka.blk; stored bond number as a 356* function of atom number and connection number 357* iob: local integer; atom number of second atom in bond 358* ncubo: local integer; stores number of current bond 359* na: stored in cbka.blk: number of atoms in system 360* zero: stored in cbka.blk: value 0.00 361* 362********************************************************************** 363 do 10 i1=1,nbon 364 ibt=ib(i1,1) 365 j1=ib(i1,2) 366 j2=ib(i1,3) 367 if (ovc(ibt).lt.0.001.and.v13cor(ibt).lt.0.001) then 368 idbo1(i1)=2 369 idbo(i1,1)=j1 370 idbo(i1,2)=j2 371 do k1=1,3 372 dbondc(i1,k1,1)=dbodc(i1,k1,1) 373 dbondc(i1,k1,2)=dbodc(i1,k1,2) 374 dbosindc(i1,k1,1)=dbosidc(i1,k1,1) 375 dbosindc(i1,k1,2)=dbosidc(i1,k1,2) 376 dbopindc(i1,k1,1)=dbopidc(i1,k1,1) 377 dbopindc(i1,k1,2)=dbopidc(i1,k1,2) 378 dbopi2ndc(i1,k1,1)=dbopi2dc(i1,k1,1) 379 dbopi2ndc(i1,k1,2)=dbopi2dc(i1,k1,2) 380 end do 381 goto 10 382 end if 383 boo=bo(i1) 384 bopio=bopi(i1) 385 bopi2o=bopi2(i1) 386 iti=ia(j1,1) 387 itj=ia(j2,1) 388 aboi=abo(j1) 389 aboj=abo(j2) 390 vp131=sqrt(bo131(iti)*bo131(itj)) 391 vp132=sqrt(bo132(iti)*bo132(itj)) 392 vp133=sqrt(bo133(iti)*bo133(itj)) 393 corrtot=1.0 394 dbodsboi1=zero 395 dbodsboj1=zero 396 if (ovc(ibt).gt.0.001) then 397 ovi=aboi-aval(iti) 398 ovj=aboj-aval(itj) 399 400********************************************************************** 401* * 402* Correction for overcoordination * 403* * 404********************************************************************** 405 exphu1=exp(-vpar(2)*ovi) 406 exphu2=exp(-vpar(2)*ovj) 407 exp11=exp(-vpar(1)*ovi) 408 exp21=exp(-vpar(1)*ovj) 409 exphu12=(exphu1+exphu2) 410 ovcor=-(1.0/vpar(2))*log(0.50*exphu12) 411* huli=((1.0/ovc(ibt))*aval(iti)+exp11+exp21) 412* hulj=((1.0/ovc(ibt))*aval(itj)+exp11+exp21) 413 huli=aval(iti)+exp11+exp21 414 hulj=aval(itj)+exp11+exp21 415 corr1=huli/(huli+ovcor) 416 corr2=hulj/(hulj+ovcor) 417 corrtot=0.50*(corr1+corr2) 418 419 dbodsboi1=0.50*(-vpar(1)*exp11/(huli+ovcor)- 420 $(corr1/(huli+ovcor))* 421 $(-vpar(1)*exp11+exphu1/exphu12)-vpar(1)*exp11/(hulj+ovcor)- 422 $(corr2/(hulj+ovcor))*(-vpar(1)*exp11+exphu1/exphu12)) 423 dbodsboj1=0.50*(-vpar(1)*exp21/(huli+ovcor)- 424 $(corr1/(huli+ovcor))* 425 $(-vpar(1)*exp21+exphu2/exphu12)-vpar(1)*exp21/(hulj+ovcor)- 426 $(corr2/(hulj+ovcor))*(-vpar(1)*exp21+exphu2/exphu12)) 427 end if 428********************************************************************** 429* * 430* Correction for 1-3 bond orders * 431* * 432********************************************************************** 433 dbodsboi2=zero 434 dbodsboj2=zero 435 bocor1=1.0 436 bocor2=1.0 437 if (v13cor(ibt).gt.0.001) then 438 ovi2=aboi-vval3(iti) !Modification for metal surfaces 439 ovj2=aboj-vval3(itj) 440* ovi2=aboi-valf(iti) 441* ovj2=aboj-valf(itj) 442* ovi2=aboi-aval(iti) 443* ovj2=aboj-aval(itj) 444 cor1=vp131*boo*boo-ovi2 445 cor2=vp131*boo*boo-ovj2 446* exphu3=v13cor(ibt)*exp(-vp132*cor1+vp133) 447* exphu4=v13cor(ibt)*exp(-vp132*cor2+vp133) 448 exphu3=exp(-vp132*cor1+vp133) 449 exphu4=exp(-vp132*cor2+vp133) 450 bocor1=1.0/(1.0+exphu3) 451 bocor2=1.0/(1.0+exphu4) 452 dbodsboi2=-bocor1*bocor1*bocor2*vp132*exphu3 453 dbodsboj2=-bocor1*bocor2*bocor2*vp132*exphu4 454 end if 455 456 bo(i1)=boo*corrtot*bocor1*bocor2 457 if (bo(i1).lt.1e-10) bo(i1)=zero 458 corrtot2=corrtot*corrtot 459 bopi(i1)=bopio*corrtot2*bocor1*bocor2 460 bopi2(i1)=bopi2o*corrtot2*bocor1*bocor2 461 if (bopi(i1).lt.1e-10) bopi(i1)=zero 462 if (bopi2(i1).lt.1e-10) bopi2(i1)=zero 463 464 dbodboo=corrtot*bocor1*bocor2+corrtot* 465 $bocor1*bocor1*bocor2*boo*vp132*vp131*2.0*boo*exphu3+ 466 $corrtot*bocor1*bocor2*bocor2*boo* 467 $vp132*vp131*exphu4*2.0*boo 468 469 dbopidbopio=corrtot2*bocor1*bocor2 470 471 dbopidboo=corrtot2* 472 $bocor1*bocor1*bocor2*boo*vp132*vp131*2.0*bopio*exphu3+ 473 $corrtot2*bocor1*bocor2*bocor2*boo* 474 $vp132*vp131*exphu4*2.0*bopio 475 476 dbopi2dbopi2o=corrtot2*bocor1*bocor2 477 478 dbopi2dboo=corrtot2* 479 $bocor1*bocor1*bocor2*boo*vp132*vp131*2.0*bopi2o*exphu3+ 480 $corrtot2*bocor1*bocor2*bocor2*boo* 481 $vp132*vp131*exphu4*2.0*bopi2o 482 483 dbodsboit=boo*dbodsboi1*bocor1*bocor2+boo*corrtot*dbodsboi2 484 dbodsbojt=boo*dbodsboj1*bocor1*bocor2+boo*corrtot*dbodsboj2 485 486 vhui=2.0*corrtot*dbodsboi1*bocor1*bocor2+corrtot2*dbodsboi2 487 vhuj=2.0*corrtot*dbodsboj1*bocor1*bocor2+corrtot2*dbodsboj2 488 dbopidsboit=bopio*vhui 489 dbopidsbojt=bopio*vhuj 490 491 dbopi2dsboit=bopi2o*vhui 492 dbopi2dsbojt=bopi2o*vhuj 493 494********************************************************************** 495* * 496* Calculate bond order derivatives * 497* * 498********************************************************************** 499 idbo1(i1)=2+ia(j1,2)+ia(j2,2) 500 idbo(i1,1)=j1 501 idbo(i1,2)=j2 502 nco=0 503 do k1=1,3 504 dbondc(i1,k1,1)=dbodc(i1,k1,1)*dbodboo 505 dbondc(i1,k1,2)=dbodc(i1,k1,2)*dbodboo 506* dbosindc(i1,k1,1)=dbosidc(i1,k1,1)*dbosidboo 507* dbosindc(i1,k1,2)=dbosidc(i1,k1,2)*dbosidboo 508 dbopindc(i1,k1,1)=dbopidc(i1,k1,1)*dbopidbopio+ 509 $dbodc(i1,k1,1)*dbopidboo 510 dbopindc(i1,k1,2)=dbopidc(i1,k1,2)*dbopidbopio+ 511 $dbodc(i1,k1,2)*dbopidboo 512 dbopi2ndc(i1,k1,1)=dbopi2dc(i1,k1,1)*dbopi2dbopi2o+ 513 $dbodc(i1,k1,1)*dbopi2dboo 514 dbopi2ndc(i1,k1,2)=dbopi2dc(i1,k1,2)*dbopi2dbopi2o+ 515 $dbodc(i1,k1,2)*dbopi2dboo 516 end do 517 do i2=1,ia(j1,2) 518 ihl=0 519 iob=ia(j1,2+i2) 520 if (iob.lt.j1) ihl=1 521 ncubo=nubon2(j1,i2) 522 idbo(i1,2+nco+1)=iob 523 do k1=1,3 524 dbondc(i1,k1,1)=dbondc(i1,k1,1)+dbodc(ncubo,k1,1+ihl)*dbodsboit 525 dbondc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbodsboit 526 527* dbosindc(i1,k1,1)=dbosindc(i1,k1,1)+ 528* $dbodc(ncubo,k1,1+ihl)*dbosidsboit 529* dbosindc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbosidsboit 530 531 dbopindc(i1,k1,1)=dbopindc(i1,k1,1)+ 532 $dbodc(ncubo,k1,1+ihl)*dbopidsboit 533 dbopindc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbopidsboit 534 535 dbopi2ndc(i1,k1,1)=dbopi2ndc(i1,k1,1)+ 536 $dbodc(ncubo,k1,1+ihl)*dbopi2dsboit 537 dbopi2ndc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbopi2dsboit 538 539 end do 540 nco=nco+1 541 end do 542 do i2=1,ia(j2,2) 543 ihl=0 544 iob=ia(j2,2+i2) 545 if (iob.lt.j2) ihl=1 546 ncubo=nubon2(j2,i2) 547 idbo(i1,2+nco+1)=iob 548 do k1=1,3 549 550 dbondc(i1,k1,2)=dbondc(i1,k1,2)+dbodc(ncubo,k1,1+ihl)*dbodsbojt 551 dbondc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbodsbojt 552 553* dbosindc(i1,k1,2)=dbosindc(i1,k1,2)+ 554* $dbodc(ncubo,k1,1+ihl)*dbosidsbojt 555* dbosindc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbosidsbojt 556 557 dbopindc(i1,k1,2)=dbopindc(i1,k1,2)+ 558 $dbodc(ncubo,k1,1+ihl)*dbopidsbojt 559 dbopindc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbopidsbojt 560 561 dbopi2ndc(i1,k1,2)=dbopi2ndc(i1,k1,2)+ 562 $dbodc(ncubo,k1,1+ihl)*dbopi2dsbojt 563 dbopi2ndc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbopi2dsbojt 564 565 end do 566 nco=nco+1 567 end do 568 569 10 continue 570 571 do i1=1,na 572 abo(i1)=zero 573 end do 574* do i1=1,na 575* do i2=1,ia(i1,2) 576* iob=ia(i1,2+i2) 577* ncubo=nubon2(i1,i2) 578* abo(i1)=abo(i1)+bo(ncubo) 579* end do 580* end do 581 do i1=1,nbon 582 j1=ib(i1,2) 583 j2=ib(i1,3) 584 abo(j1)=abo(j1)+bo(i1) 585 if (j1.ne.j2) abo(j2)=abo(j2)+bo(i1) 586 end do 587 588 15 continue 589 return 590 end 591********************************************************************** 592********************************************************************** 593 594 subroutine lonpar 595 596********************************************************************** 597#include "cbka.blk" 598#include "cbkabo.blk" 599#include "cbkconst.blk" 600#include "cbkc.blk" 601#include "cbkd.blk" 602#include "cbkdcell.blk" 603#include "cbkenergies.blk" 604#include "cbkff.blk" 605#include "cbkia.blk" 606#include "cbkidbo.blk" 607#include "cbklonpar.blk" 608#include "cbknubon2.blk" 609#include "control.blk" 610#include "small.blk" 611 dimension virial_tmp(3,3),virialsym(6) 612c$$$ if (ndebug.eq.1) then 613c$$$C open (65,file='fort.65',status='unknown',access='append') 614c$$$ write (65,*) 'In lonpar' 615c$$$ call timer(65) 616c$$$ close (65) 617c$$$ end if 618********************************************************************** 619* * 620* Calculate lone pair energy and first derivatives * 621* * 622********************************************************************** 623 elp=zero 624 do i1=1,na 625********************************************************************** 626* * 627* Determine number of lone pairs on atoms 628* * 629********************************************************************** 630 ity=ia(i1,1) 631 voptlp=0.50*(stlp(ity)-aval(ity)) 632 vlp(i1)=zero 633 vund=abo(i1)-stlp(ity) 634 vlph=2.0*int(vund/2.0) 635 vlpex=vund-vlph 636 vp16h=vpar(16)-1.0 637 638 expvlp=exp(-vpar(16)*(2.0+vlpex)*(2.0+vlpex)) 639 dvlpdsbo(i1)=-vpar(16)*2.0*(2.0+vlpex)*expvlp 640 vlp(i1)=expvlp-int(vund/2.0) 641* expvlp=exp(-vpar(16)*(2.0+vlpex)) 642* dvlpdsbo(i1)=-vpar(16)*expvlp 643* expvlp=exp(-6.0*((-0.50*vlpex)**vpar(16))) 644* vlp(i1)=(1.0-expvlp)-int(vund/2.0) 645* dvlpdsbo(i1)=-0.5*6.0*vpar(16)*((-0.5*vlpex)**vp16h)* 646* $expvlp 647********************************************************************** 648* * 649* Calculate lone pair energy * 650* * 651********************************************************************** 652 if (i1 .le. na_local) then 653 654 diffvlp=voptlp-vlp(i1) 655 exphu1=exp(-75.0*diffvlp) 656 hulp1=1.0/(1.0+exphu1) 657 elph=vlp1(ity)*diffvlp*hulp1 658* elph=vlp1(ity)*diffvlp 659 delpdvlp=-vlp1(ity)*hulp1-vlp1(ity)*diffvlp*hulp1*hulp1* 660 $75.0*exphu1 661 662 elp=elp+elph 663 estrain(i1)=estrain(i1)+elph !atom energy 664 665 delpdsbo=delpdvlp*dvlpdsbo(i1) 666********************************************************************** 667* * 668* Calculate first derivative of lone pair energy to * 669* cartesian coordinates * 670* * 671********************************************************************** 672 do i3=1,ia(i1,2) 673 iob=ia(i1,2+i3) 674 ncubo=nubon2(i1,i3) 675 676 if (Lvirial.eq.1) then 677 do k1=1,3 678 do k2=1,3 679 virial_tmp(k1,k2) = 0.0 680 end do 681 end do 682 endif 683 684 do i4=1,idbo1(ncubo) 685 ihu=idbo(ncubo,i4) 686 do k1=1,3 687 ftmp = delpdsbo*dbondc(ncubo,k1,i4) 688 d(k1,ihu)=d(k1,ihu)+ftmp 689 690 if (Lvirial.eq.1) then 691 do k1p=1,3 692 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 693 end do 694 endif 695 696 end do 697 end do 698 if (Lvirial.eq.1) then 699 virialsym(1) = virial_tmp(1,1) 700 virialsym(2) = virial_tmp(2,2) 701 virialsym(3) = virial_tmp(3,3) 702 virialsym(4) = virial_tmp(1,2) 703 virialsym(5) = virial_tmp(1,3) 704 virialsym(6) = virial_tmp(2,3) 705 do k1 = 1,6 706 virial(k1) = virial(k1) + virialsym(k1) 707 end do 708 709 if (Latomvirial.eq.1) then 710 frac = 1.0d0/idbo1(ncubo) 711 do k1 = 1,6 712 vtmp = virialsym(k1)*frac 713 do i4=1,idbo1(ncubo) 714 ihu=idbo(ncubo,i4) 715 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 716 end do 717 end do 718 endif 719 endif 720 721 end do 722 endif 723 724 end do 725 726 return 727 end 728********************************************************************** 729********************************************************************** 730 731 subroutine covbon 732 733********************************************************************** 734#include "cbka.blk" 735#include "cbkc.blk" 736#include "cbkabo.blk" 737#include "cbkbo.blk" 738#include "cbkbosi.blk" 739#include "cbkbopi.blk" 740#include "cbkbopi2.blk" 741#include "cbkconst.blk" 742#include "cbkcovbon.blk" 743#include "cbkd.blk" 744#include "cbkdbopi2ndc.blk" 745#include "cbkdbopindc.blk" 746#include "cbkdcell.blk" 747#include "cbkenergies.blk" 748#include "cbkff.blk" 749#include "cbkia.blk" 750#include "cbkidbo.blk" 751#include "cbknubon2.blk" 752#include "cbkqa.blk" 753#include "cbkrbo.blk" 754#include "control.blk" 755#include "small.blk" 756 dimension virial_tmp(3,3),virialsym(6) 757********************************************************************** 758* * 759* Calculate bond energy and first derivatives * 760* * 761********************************************************************** 762c$$$ if (ndebug.eq.1) then 763c$$$C open (65,file='fort.65',status='unknown',access='append') 764c$$$ write (65,*) 'In covbon' 765c$$$ call timer(65) 766c$$$ close (65) 767c$$$ end if 768 eb=0.0d0 769 if (nbon.eq.0) return 770********************************************************************** 771* * 772* Calculate bond energies * 773* * 774********************************************************************** 775c$$$ if (ndebug.eq.1) then 776c$$$C open (65,file='fort.65',status='unknown',access='append') 777c$$$ write(65,*) 'Bond forces' 778c$$$ write(65,*) 'nbon = ',nbon 779c$$$ endif 780 781 do 20 i1=1,nbon 782 783 boa=bo(i1) 784* if (boa.lt.cutof2) goto 20 785 j1=ib(i1,2) 786 j2=ib(i1,3) 787 788c Only compute interaction if both atoms 789c are local or else flip a coin 790 if (j1 .gt. na_local) go to 20 791 if (j2 .gt. na_local) then 792 if (itag(j1) .lt. itag(j2)) go to 20 793 if (itag(j1) .eq. itag(j2)) then 794 if(c(j1,3) .gt. c(j2,3)) go to 20 795 if(c(j1,3) .eq. c(j2,3) .and. 796 $ c(j1,2) .gt. c(j2,2)) go to 20 797 if(c(j1,3) .eq. c(j2,3) .and. 798 $ c(j1,2) .eq. c(j2,2) .and. 799 $ c(j1,1) .gt. c(j2,1)) go to 20 800 endif 801 endif 802 vsymm=1.0 803 if (j1.eq.j2) vsymm=0.5 804 805 bopia=bopi(i1) 806 bopi2a=bopi2(i1) 807 bosia=boa-bopia-bopi2a 808 if (bosia.lt.zero) bosia=zero 809 it1=ia(j1,1) 810 it2=ia(j2,1) 811 ibt=ib(i1,1) 812 de1h=vsymm*de1(ibt) 813 de2h=vsymm*de2(ibt) 814 de3h=vsymm*de3(ibt) 815 816 bopo1=bosia**psp(ibt) 817 exphu1=exp(psi(ibt)*(1.0-bopo1)) 818 ebh=-de1h*bosia*exphu1-de2h*bopia-de3h*bopi2a 819 820 debdbo=-de1h*exphu1+de1h*exphu1*psp(ibt)*psi(ibt)*bopo1 821 debdbopi=-de2h 822 debdbopi2=-de3h 823 824 eb=eb+ebh 825 estrain(j1)=estrain(j1)+0.50*ebh !1st atom energy 826 estrain(j2)=estrain(j2)+0.50*ebh !2nd atom energy 827 828 if (Lvirial.eq.1) then 829 do k1=1,3 830 do k2=1,3 831 virial_tmp(k1,k2) = 0.0 832 end do 833 end do 834 endif 835 836 do i2=1,idbo1(i1) 837 ihu=idbo(i1,i2) 838 do k1=1,3 839 ftmp = debdbo*(dbondc(i1,k1,i2)-dbopindc(i1,k1,i2)- 840 $dbopi2ndc(i1,k1,i2))+ 841 $debdbopi*dbopindc(i1,k1,i2)+ 842 $debdbopi2*dbopi2ndc(i1,k1,i2) 843 d(k1,ihu)=d(k1,ihu)+ftmp 844 845 if (Lvirial.eq.1) then 846 do k1p=1,3 847 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 848 end do 849 endif 850 851 end do 852 end do 853 if (Lvirial.eq.1) then 854 virialsym(1) = virial_tmp(1,1) 855 virialsym(2) = virial_tmp(2,2) 856 virialsym(3) = virial_tmp(3,3) 857 virialsym(4) = virial_tmp(1,2) 858 virialsym(5) = virial_tmp(1,3) 859 virialsym(6) = virial_tmp(2,3) 860 do k1 = 1,6 861 virial(k1) = virial(k1) + virialsym(k1) 862 end do 863 864 if (Latomvirial.eq.1) then 865 frac = 1.0d0/idbo1(i1) 866 do k1 = 1,6 867 vtmp = virialsym(k1)*frac 868 do i2=1,idbo1(i1) 869 ihu=idbo(i1,i2) 870 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 871 end do 872 end do 873 endif 874 875 endif 876 877********************************************************************** 878* * 879* Stabilisation terminal triple bond in CO * 880* * 881********************************************************************** 882 if (boa.lt.1.00) goto 20 883* Stabilization for all triple bonds (not just for CO) in ReaxFF combustion FF 884 if (ltripstaball.eq.1 .or. 885 $ (qa(j1).eq.'C '.and.qa(j2).eq.'O ').or. 886 $ (qa(j1).eq.'O '.and.qa(j2).eq.'C ')) then 887 888 ba=(boa-2.50)*(boa-2.50) 889 exphu=exp(-vpar(8)*ba) 890 oboa=abo(j1)-boa 891 obob=abo(j2)-boa 892 exphua1=exp(-vpar(4)*oboa) 893 exphub1=exp(-vpar(4)*obob) 894 ovoab=abo(j1)-aval(it1)+abo(j2)-aval(it2) 895 exphuov=exp(vpar(5)*ovoab) 896 hulpov=1.0/(1.0+25.0*exphuov) 897 898 estriph=vpar(11)*exphu*hulpov*(exphua1+exphub1) 899 900 eb=eb+estriph 901 estrain(j1)=estrain(j1)+0.50*estriph !1st atom energy 902 estrain(j2)=estrain(j2)+0.50*estriph !2nd atom energy 903 904 decobdbo=vpar(4)*vpar(11)*exphu*hulpov*(exphua1+exphub1) 905 $-2.0*vpar(11)*vpar(8)*(boa-2.50)*hulpov*exphu* 906 $(exphua1+exphub1) 907 decobdboua=-25.0*vpar(5)*vpar(11)*exphu*exphuov*hulpov*hulpov* 908 $(exphua1+exphub1)-vpar(11)*exphu*vpar(4)*hulpov*exphua1 909 decobdboub=-25.0*vpar(5)*vpar(11)*exphu*exphuov*hulpov*hulpov* 910 $(exphua1+exphub1)-vpar(11)*exphu*vpar(4)*hulpov*exphub1 911 912 if (Lvirial.eq.1) then 913 do k1=1,3 914 do k2=1,3 915 virial_tmp(k1,k2) = 0.0 916 end do 917 end do 918 endif 919 920 do i2=1,idbo1(i1) 921 ihu=idbo(i1,i2) 922 do k1=1,3 923 ftmp = decobdbo*dbondc(i1,k1,i2) 924 d(k1,ihu)=d(k1,ihu)+ftmp 925 926 if (Lvirial.eq.1) then 927 do k1p=1,3 928 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 929 end do 930 endif 931 932 end do 933 end do 934 if (Lvirial.eq.1) then 935 virialsym(1) = virial_tmp(1,1) 936 virialsym(2) = virial_tmp(2,2) 937 virialsym(3) = virial_tmp(3,3) 938 virialsym(4) = virial_tmp(1,2) 939 virialsym(5) = virial_tmp(1,3) 940 virialsym(6) = virial_tmp(2,3) 941 do k1 = 1,6 942 virial(k1) = virial(k1) + virialsym(k1) 943 end do 944 945 if (Latomvirial.eq.1) then 946 frac = 1.0d0/idbo1(i1) 947 do k1 = 1,6 948 vtmp = virialsym(k1)*frac 949 do i2=1,idbo1(i1) 950 ihu=idbo(i1,i2) 951 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 952 end do 953 end do 954 endif 955 956 endif 957 958 do i3=1,ia(j1,2) 959 iob=ia(j1,2+i3) 960 ncubo=nubon2(j1,i3) 961 962 if (Lvirial.eq.1) then 963 do k1=1,3 964 do k2=1,3 965 virial_tmp(k1,k2) = 0.0 966 end do 967 end do 968 endif 969 970 do i4=1,idbo1(ncubo) 971 ihu=idbo(ncubo,i4) 972 do k1=1,3 973 ftmp = decobdboua*dbondc(ncubo,k1,i4) 974 d(k1,ihu)=d(k1,ihu)+ftmp 975 976 if (Lvirial.eq.1) then 977 do k1p=1,3 978 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 979 end do 980 endif 981 982 end do 983 end do 984 if (Lvirial.eq.1) then 985 virialsym(1) = virial_tmp(1,1) 986 virialsym(2) = virial_tmp(2,2) 987 virialsym(3) = virial_tmp(3,3) 988 virialsym(4) = virial_tmp(1,2) 989 virialsym(5) = virial_tmp(1,3) 990 virialsym(6) = virial_tmp(2,3) 991 do k1 = 1,6 992 virial(k1) = virial(k1) + virialsym(k1) 993 end do 994 995 if (Latomvirial.eq.1) then 996 frac = 1.0d0/idbo1(ncubo) 997 do k1 = 1,6 998 vtmp = virialsym(k1)*frac 999 do i4=1,idbo1(ncubo) 1000 ihu=idbo(ncubo,i4) 1001 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 1002 end do 1003 end do 1004 endif 1005 1006 endif 1007 1008 end do 1009 1010 do i3=1,ia(j2,2) 1011 iob=ia(j2,2+i3) 1012 ncubo=nubon2(j2,i3) 1013 if (Lvirial.eq.1) then 1014 do k1=1,3 1015 do k2=1,3 1016 virial_tmp(k1,k2) = 0.0 1017 end do 1018 end do 1019 endif 1020 do i4=1,idbo1(ncubo) 1021 ihu=idbo(ncubo,i4) 1022 do k1=1,3 1023 ftmp = decobdboub*dbondc(ncubo,k1,i4) 1024 d(k1,ihu)=d(k1,ihu)+ftmp 1025 1026 if (Lvirial.eq.1) then 1027 do k1p=1,3 1028 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 1029 end do 1030 endif 1031 1032 end do 1033 end do 1034 if (Lvirial.eq.1) then 1035 virialsym(1) = virial_tmp(1,1) 1036 virialsym(2) = virial_tmp(2,2) 1037 virialsym(3) = virial_tmp(3,3) 1038 virialsym(4) = virial_tmp(1,2) 1039 virialsym(5) = virial_tmp(1,3) 1040 virialsym(6) = virial_tmp(2,3) 1041 do k1 = 1,6 1042 virial(k1) = virial(k1) + virialsym(k1) 1043 end do 1044 1045 if (Latomvirial.eq.1) then 1046 frac = 1.0d0/idbo1(ncubo) 1047 do k1 = 1,6 1048 vtmp = virialsym(k1)*frac 1049 do i4=1,idbo1(ncubo) 1050 ihu=idbo(ncubo,i4) 1051 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 1052 end do 1053 end do 1054 endif 1055 1056 endif 1057 1058 end do 1059 1060 endif 1061 1062 20 continue 1063 1064 return 1065 end 1066********************************************************************** 1067********************************************************************** 1068 1069 subroutine ovcor 1070 1071********************************************************************** 1072#include "cbka.blk" 1073#include "cbkc.blk" 1074#include "cbkabo.blk" 1075#include "cbkbo.blk" 1076#include "cbkbopi.blk" 1077#include "cbkbopi2.blk" 1078#include "cbkconst.blk" 1079#include "cbkd.blk" 1080#include "cbkdbopi2ndc.blk" 1081#include "cbkdbopindc.blk" 1082#include "cbkdcell.blk" 1083#include "cbkenergies.blk" 1084#include "cbkff.blk" 1085#include "cbkia.blk" 1086#include "cbkidbo.blk" 1087#include "cbklonpar.blk" 1088#include "cbknubon2.blk" 1089#include "cbkrbo.blk" 1090#include "control.blk" 1091#include "small.blk" 1092********************************************************************** 1093* * 1094* Calculate atom energy * 1095* Correction for over- and undercoordinated atoms * 1096* * 1097********************************************************************** 1098 dimension vlptemp(nat) 1099 dimension virial_tmp(3,3),virialsym(6) 1100c$$$ if (ndebug.eq.1) then 1101c$$$C open (65,file='fort.65',status='unknown',access='append') 1102c$$$ write (65,*) 'In ovcor' 1103c$$$ call timer(65) 1104c$$$ close (65) 1105c$$$ end if 1106 do i1=1,na 1107 ity1=ia(i1,1) 1108 vlptemp(i1)=vlp(i1) 1109 if (amas(ity1).gt.21.0) vlptemp(i1)=0.50*(stlp(ity1)-aval(ity1)) !Only for 1st-row elements 1110 end do 1111 25 ea=zero 1112 eaot=zero 1113 eaut=zero 1114 epen=0.0 1115 1116 do 30 i1=1,na_local 1117 ity1=ia(i1,1) 1118 dfvl=1.0 1119 if (amas(ity1).gt.21.0) dfvl=0.0 !Only for 1st-row elements 1120********************************************************************** 1121* * 1122* Calculate overcoordination energy * 1123* Valency is corrected for lone pairs * 1124* * 1125********************************************************************** 1126 1127 voptlp=0.50*(stlp(ity1)-aval(ity1)) 1128 diffvlph=dfvl*(voptlp-vlptemp(i1)) 1129********************************************************************** 1130* * 1131* Determine coordination neighboring atoms * 1132* * 1133********************************************************************** 1134 sumov=0.0 1135 sumov2=0.0 1136 do i3=1,ia(i1,2) 1137 iat2=ia(i1,2+i3) 1138 ity2=ia(iat2,1) 1139 ncubo=nubon2(i1,i3) 1140 if (bo(ncubo).gt.0.0) then 1141 ibt=ib(ncubo,1) 1142 voptlp2=0.50*(stlp(ity2)-aval(ity2)) 1143 diffvlp2=dfvl*(voptlp2-vlptemp(iat2)) 1144 sumov=sumov+(bopi(ncubo)+bopi2(ncubo))* 1145 $(abo(iat2)-aval(ity2)-diffvlp2) 1146 sumov2=sumov2+vover(ibt)*de1(ibt)*bo(ncubo) 1147 endif 1148 end do 1149 1150 exphu1=exp(vpar(32)*sumov) 1151 vho=1.0/(1.0+vpar(33)*exphu1) 1152 diffvlp=diffvlph*vho 1153 1154 vov1=abo(i1)-aval(ity1)-diffvlp 1155 dvov1dsumov=diffvlph*vpar(32)*vpar(33)*vho*vho*exphu1 1156 exphuo=exp(vovun(ity1)*vov1) 1157 hulpo=1.0/(1.0+exphuo) 1158 1159 hulpp=(1.0/(vov1+aval(ity1)+1e-8)) 1160 1161 eah=sumov2*hulpp*hulpo*vov1 1162 deadvov1=-sumov2*hulpp*hulpp*vov1*hulpo+ 1163 $sumov2*hulpp*hulpo-sumov2*hulpp*vov1*vovun(ity1)* 1164 $hulpo*hulpo*exphuo 1165 1166 ea=ea+eah 1167 estrain(i1)=estrain(i1)+eah !atom energy 1168********************************************************************** 1169* * 1170* Calculate first derivative of overcoordination energy to * 1171* cartesian coordinates * 1172* * 1173********************************************************************** 1174 do i3=1,ia(i1,2) 1175 iob=ia(i1,2+i3) 1176 ncubo=nubon2(i1,i3) 1177 if (bo(ncubo).gt.0.0) then 1178 ibt=ib(ncubo,1) 1179 deadbo=vover(ibt)*de1(ibt)*hulpp*hulpo*vov1 1180 1181 if (Lvirial.eq.1) then 1182 do k1=1,3 1183 do k2=1,3 1184 virial_tmp(k1,k2) = 0.0 1185 end do 1186 end do 1187 endif 1188 1189 do i4=1,idbo1(ncubo) 1190 ihu=idbo(ncubo,i4) 1191 do k1=1,3 1192 ftmp = deadvov1*(1.0+dfvl*vho*dvlpdsbo(i1))* 1193 $dbondc(ncubo,k1,i4)+deadbo*dbondc(ncubo,k1,i4) 1194 d(k1,ihu)=d(k1,ihu)+ftmp 1195 1196 if (Lvirial.eq.1) then 1197 do k1p=1,3 1198 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 1199 end do 1200 endif 1201 1202 end do 1203 end do 1204 if (Lvirial.eq.1) then 1205 virialsym(1) = virial_tmp(1,1) 1206 virialsym(2) = virial_tmp(2,2) 1207 virialsym(3) = virial_tmp(3,3) 1208 virialsym(4) = virial_tmp(1,2) 1209 virialsym(5) = virial_tmp(1,3) 1210 virialsym(6) = virial_tmp(2,3) 1211 do k1 = 1,6 1212 virial(k1) = virial(k1) + virialsym(k1) 1213 end do 1214 1215 if (Latomvirial.eq.1) then 1216 frac = 1.0d0/idbo1(ncubo) 1217 do k1 = 1,6 1218 vtmp = virialsym(k1)*frac 1219 do i4=1,idbo1(ncubo) 1220 ihu=idbo(ncubo,i4) 1221 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 1222 end do 1223 end do 1224 endif 1225 1226 endif 1227 1228 endif 1229 end do 1230 1231 do i2=1,ia(i1,2) 1232 1233 iat2=ia(i1,2+i2) 1234 ity2=ia(iat2,1) 1235 nbosa=nubon2(i1,i2) 1236 if (bo(nbosa).gt.0.0) then 1237 deadvov2=deadvov1*dvov1dsumov*(bopi(nbosa)+bopi2(nbosa)) 1238 1239 voptlp2=0.50*(stlp(ity2)-aval(ity2)) 1240 diffvlp2=dfvl*(voptlp2-vlptemp(iat2)) 1241 deadpibo=deadvov1*dvov1dsumov*(abo(iat2)-aval(ity2)-diffvlp2) 1242 1243 if (Lvirial.eq.1) then 1244 do k1=1,3 1245 do k2=1,3 1246 virial_tmp(k1,k2) = 0.0 1247 end do 1248 end do 1249 endif 1250 1251 do i4=1,idbo1(nbosa) 1252 ihu=idbo(nbosa,i4) 1253 do k1=1,3 1254 ftmp = deadpibo*(dbopindc(nbosa,k1,i4)+ 1255 $dbopi2ndc(nbosa,k1,i4)) 1256 d(k1,ihu)=d(k1,ihu)+ftmp 1257 1258 if (Lvirial.eq.1) then 1259 do k1p=1,3 1260 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 1261 end do 1262 endif 1263 1264 end do 1265 end do 1266 if (Lvirial.eq.1) then 1267 virialsym(1) = virial_tmp(1,1) 1268 virialsym(2) = virial_tmp(2,2) 1269 virialsym(3) = virial_tmp(3,3) 1270 virialsym(4) = virial_tmp(1,2) 1271 virialsym(5) = virial_tmp(1,3) 1272 virialsym(6) = virial_tmp(2,3) 1273 do k1 = 1,6 1274 virial(k1) = virial(k1) + virialsym(k1) 1275 end do 1276 1277 if (Latomvirial.eq.1) then 1278 frac = 1.0d0/idbo1(nbosa) 1279 do k1 = 1,6 1280 vtmp = virialsym(k1)*frac 1281 do i4=1,idbo1(nbosa) 1282 ihu=idbo(nbosa,i4) 1283 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 1284 end do 1285 end do 1286 endif 1287 endif 1288 1289 do i3=1,ia(iat2,2) 1290 iob=ia(iat2,2+i3) 1291 ncubo=nubon2(iat2,i3) 1292 if (bo(ncubo).gt.0.0) then 1293 1294 if (Lvirial.eq.1) then 1295 do k1=1,3 1296 do k2=1,3 1297 virial_tmp(k1,k2) = 0.0 1298 end do 1299 end do 1300 endif 1301 1302 do i4=1,idbo1(ncubo) 1303 ihu=idbo(ncubo,i4) 1304 do k1=1,3 1305 ftmp = deadvov2*(1.0+dfvl*dvlpdsbo(iat2))* 1306 $dbondc(ncubo,k1,i4) 1307 d(k1,ihu)=d(k1,ihu)+ftmp 1308 1309 if (Lvirial.eq.1) then 1310 do k1p=1,3 1311 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 1312 end do 1313 endif 1314 1315 end do 1316 end do 1317 if (Lvirial.eq.1) then 1318 virialsym(1) = virial_tmp(1,1) 1319 virialsym(2) = virial_tmp(2,2) 1320 virialsym(3) = virial_tmp(3,3) 1321 virialsym(4) = virial_tmp(1,2) 1322 virialsym(5) = virial_tmp(1,3) 1323 virialsym(6) = virial_tmp(2,3) 1324 do k1 = 1,6 1325 virial(k1) = virial(k1) + virialsym(k1) 1326 end do 1327 1328 if (Latomvirial.eq.1) then 1329 frac = 1.0d0/idbo1(ncubo) 1330 do k1 = 1,6 1331 vtmp = virialsym(k1)*frac 1332 do i4=1,idbo1(ncubo) 1333 ihu=idbo(ncubo,i4) 1334 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 1335 end do 1336 end do 1337 endif 1338 1339 endif 1340 1341 endif 1342 end do 1343 1344 endif 1345 1346 end do 1347 1348********************************************************************** 1349* * 1350* Calculate undercoordination energy * 1351* * 1352********************************************************************** 1353 if (valp1(ity1).lt.zero) goto 30 !skip undercoordination 1354 exphu2=exp(vpar(10)*sumov) 1355 vuhu1=1.0+vpar(9)*exphu2 1356 hulpu2=1.0/vuhu1 1357 1358 exphu3=-exp(vpar(7)*vov1) 1359 hulpu3=-(1.0+exphu3) 1360 1361 dise2=valp1(ity1) 1362 exphuu=exp(-vovun(ity1)*vov1) 1363 hulpu=1.0/(1.0+exphuu) 1364 eahu=dise2*hulpu*hulpu2*hulpu3 1365 deaudvov1=dise2*hulpu2*vovun(ity1)*hulpu*hulpu*exphuu*hulpu3- 1366 $dise2*hulpu*hulpu2*vpar(7)*exphu3 1367 1368 ea=ea+eahu 1369 estrain(i1)=estrain(i1)+eahu !atom energy 1370 1371 deaudsumov=-dise2*hulpu*vpar(9)*vpar(10)*hulpu3*exphu2* 1372 $hulpu2*hulpu2 1373 1374********************************************************************** 1375* * 1376* Calculate first derivative of atom energy to cartesian * 1377* coordinates * 1378* * 1379********************************************************************** 1380 1381 do i3=1,ia(i1,2) 1382 iob=ia(i1,2+i3) 1383 ncubo=nubon2(i1,i3) 1384 if (bo(ncubo).gt.0.0) then 1385 1386 if (Lvirial.eq.1) then 1387 do k1=1,3 1388 do k2=1,3 1389 virial_tmp(k1,k2) = 0.0 1390 end do 1391 end do 1392 endif 1393 1394 do i4=1,idbo1(ncubo) 1395 ihu=idbo(ncubo,i4) 1396 do k1=1,3 1397 ftmp = deaudvov1*(1.0+dfvl*vho*dvlpdsbo(i1))* 1398 $dbondc(ncubo,k1,i4) 1399 d(k1,ihu)=d(k1,ihu)+ftmp 1400 1401 if (Lvirial.eq.1) then 1402 do k1p=1,3 1403 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 1404 end do 1405 endif 1406 1407 end do 1408 end do 1409 if (Lvirial.eq.1) then 1410 virialsym(1) = virial_tmp(1,1) 1411 virialsym(2) = virial_tmp(2,2) 1412 virialsym(3) = virial_tmp(3,3) 1413 virialsym(4) = virial_tmp(1,2) 1414 virialsym(5) = virial_tmp(1,3) 1415 virialsym(6) = virial_tmp(2,3) 1416 do k1 = 1,6 1417 virial(k1) = virial(k1) + virialsym(k1) 1418 end do 1419 1420 if (Latomvirial.eq.1) then 1421 frac = 1.0d0/idbo1(ncubo) 1422 do k1 = 1,6 1423 vtmp = virialsym(k1)*frac 1424 do i4=1,idbo1(ncubo) 1425 ihu=idbo(ncubo,i4) 1426 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 1427 end do 1428 end do 1429 endif 1430 1431 endif 1432 1433 endif 1434 end do 1435 1436 do i2=1,ia(i1,2) 1437 1438 iat2=ia(i1,2+i2) 1439 ity2=ia(iat2,1) 1440 nbosa=nubon2(i1,i2) 1441 if (bo(nbosa).gt.0.0) then 1442 deadvov2=(deaudsumov+dvov1dsumov*deaudvov1)* 1443 $(bopi(nbosa)+bopi2(nbosa)) 1444 1445 voptlp2=0.50*(stlp(ity2)-aval(ity2)) 1446 diffvlp2=dfvl*(voptlp2-vlptemp(iat2)) 1447 deadpibo1=(dvov1dsumov*deaudvov1+deaudsumov)* 1448 $(abo(iat2)-aval(ity2)-diffvlp2) 1449 1450 if (Lvirial.eq.1) then 1451 do k1=1,3 1452 do k2=1,3 1453 virial_tmp(k1,k2) = 0.0 1454 end do 1455 end do 1456 endif 1457 1458 do i4=1,idbo1(nbosa) 1459 ihu=idbo(nbosa,i4) 1460 do k1=1,3 1461 ftmp = deadpibo1* 1462 $(dbopindc(nbosa,k1,i4)+dbopi2ndc(nbosa,k1,i4)) 1463 d(k1,ihu)=d(k1,ihu)+ftmp 1464 1465 if (Lvirial.eq.1) then 1466 do k1p=1,3 1467 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 1468 end do 1469 endif 1470 1471 end do 1472 end do 1473 if (Lvirial.eq.1) then 1474 virialsym(1) = virial_tmp(1,1) 1475 virialsym(2) = virial_tmp(2,2) 1476 virialsym(3) = virial_tmp(3,3) 1477 virialsym(4) = virial_tmp(1,2) 1478 virialsym(5) = virial_tmp(1,3) 1479 virialsym(6) = virial_tmp(2,3) 1480 do k1 = 1,6 1481 virial(k1) = virial(k1) + virialsym(k1) 1482 end do 1483 1484 if (Latomvirial.eq.1) then 1485 frac = 1.0d0/idbo1(nbosa) 1486 do k1 = 1,6 1487 vtmp = virialsym(k1)*frac 1488 do i4=1,idbo1(nbosa) 1489 ihu=idbo(nbosa,i4) 1490 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 1491 end do 1492 end do 1493 endif 1494 1495 endif 1496 1497 do i3=1,ia(iat2,2) 1498 iob=ia(iat2,2+i3) 1499 ncubo=nubon2(iat2,i3) 1500 if (bo(ncubo).gt.0.0) then 1501 1502 if (Lvirial.eq.1) then 1503 do k1=1,3 1504 do k2=1,3 1505 virial_tmp(k1,k2) = 0.0 1506 end do 1507 end do 1508 endif 1509 1510 do i4=1,idbo1(ncubo) 1511 ihu=idbo(ncubo,i4) 1512 do k1=1,3 1513 ftmp = deadvov2*(1.0+dfvl*dvlpdsbo(iat2))* 1514 $dbondc(ncubo,k1,i4) 1515 d(k1,ihu)=d(k1,ihu)+ftmp 1516 1517 if (Lvirial.eq.1) then 1518 do k1p=1,3 1519 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 1520 end do 1521 endif 1522 1523 end do 1524 end do 1525 if (Lvirial.eq.1) then 1526 virialsym(1) = virial_tmp(1,1) 1527 virialsym(2) = virial_tmp(2,2) 1528 virialsym(3) = virial_tmp(3,3) 1529 virialsym(4) = virial_tmp(1,2) 1530 virialsym(5) = virial_tmp(1,3) 1531 virialsym(6) = virial_tmp(2,3) 1532 do k1 = 1,6 1533 virial(k1) = virial(k1) + virialsym(k1) 1534 end do 1535 1536 if (Latomvirial.eq.1) then 1537 frac = 1.0d0/idbo1(ncubo) 1538 do k1 = 1,6 1539 vtmp = virialsym(k1)*frac 1540 do i4=1,idbo1(ncubo) 1541 ihu=idbo(ncubo,i4) 1542 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 1543 end do 1544 end do 1545 endif 1546 1547 endif 1548 1549 endif 1550 end do 1551 1552 endif 1553 1554 end do 1555 1556 1557 30 continue 1558 1559********************************************************************** 1560* * 1561* Calculate correction for C2 * 1562* * 1563********************************************************************** 1564 if (abs(vpar(6)).gt.0.001) then 1565 do 40 i1=1,na_local 1566 ity1=ia(i1,1) 1567 vov4=abo(i1)-aval(ity1) 1568 1569 do i2=1,ia(i1,2) 1570 iat2=ia(i1,2+i2) 1571 nbohu=nubon2(i1,i2) 1572 if (bo(nbohu).gt.0.0) then 1573 1574 ibt=ib(nbohu,1) 1575 elph=zero 1576 deahu2dbo=zero 1577 deahu2dsbo=zero 1578 vov3=bo(nbohu)-vov4-0.040*(vov4**4) 1579 if (vov3.gt.3.0) then 1580 elph=vpar(6)*(vov3-3.0)*(vov3-3.0) 1581 deahu2dbo=2.0*vpar(6)*(vov3-3.0) 1582 deahu2dsbo=2.0*vpar(6)*(vov3-3.0)*(-1.0- 1583 $0.16*(vov4**3)) 1584 end if 1585 1586 elp=elp+elph 1587 estrain(i1)=estrain(i1)+elph !atom energy 1588 1589 if (Lvirial.eq.1) then 1590 do k1=1,3 1591 do k2=1,3 1592 virial_tmp(k1,k2) = 0.0 1593 end do 1594 end do 1595 endif 1596 1597 do i3=1,idbo1(nbohu) 1598 ihu=idbo(nbohu,i3) 1599 do k1=1,3 1600 ftmp = deahu2dbo*dbondc(nbohu,k1,i3) 1601 d(k1,ihu)=d(k1,ihu)+ftmp 1602 1603 if (Lvirial.eq.1) then 1604 do k1p=1,3 1605 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 1606 end do 1607 endif 1608 1609 end do 1610 end do 1611 if (Lvirial.eq.1) then 1612 virialsym(1) = virial_tmp(1,1) 1613 virialsym(2) = virial_tmp(2,2) 1614 virialsym(3) = virial_tmp(3,3) 1615 virialsym(4) = virial_tmp(1,2) 1616 virialsym(5) = virial_tmp(1,3) 1617 virialsym(6) = virial_tmp(2,3) 1618 do k1 = 1,6 1619 virial(k1) = virial(k1) + virialsym(k1) 1620 end do 1621 1622 if (Latomvirial.eq.1) then 1623 frac = 1.0d0/idbo1(nbohu) 1624 do k1 = 1,6 1625 vtmp = virialsym(k1)*frac 1626 do i3=1,idbo1(nbohu) 1627 ihu=idbo(nbohu,i3) 1628 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 1629 end do 1630 end do 1631 endif 1632 1633 endif 1634 1635 do i3=1,ia(i1,2) 1636 iob=ia(i1,2+i3) 1637 ncubo=nubon2(i1,i3) 1638 if (bo(ncubo).gt.0.0) then 1639 1640 if (Lvirial.eq.1) then 1641 do k1=1,3 1642 do k2=1,3 1643 virial_tmp(k1,k2) = 0.0 1644 end do 1645 end do 1646 endif 1647 1648 do i4=1,idbo1(ncubo) 1649 ihu=idbo(ncubo,i4) 1650 do k1=1,3 1651 ftmp = deahu2dsbo*dbondc(ncubo,k1,i4) 1652 d(k1,ihu)=d(k1,ihu)+ftmp 1653 1654 if (Lvirial.eq.1) then 1655 do k1p=1,3 1656 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 1657 end do 1658 endif 1659 1660 end do 1661 end do 1662 if (Lvirial.eq.1) then 1663 virialsym(1) = virial_tmp(1,1) 1664 virialsym(2) = virial_tmp(2,2) 1665 virialsym(3) = virial_tmp(3,3) 1666 virialsym(4) = virial_tmp(1,2) 1667 virialsym(5) = virial_tmp(1,3) 1668 virialsym(6) = virial_tmp(2,3) 1669 do k1 = 1,6 1670 virial(k1) = virial(k1) + virialsym(k1) 1671 end do 1672 1673 if (Latomvirial.eq.1) then 1674 frac = 1.0d0/idbo1(ncubo) 1675 do k1 = 1,6 1676 vtmp = virialsym(k1)*frac 1677 do i4=1,idbo1(ncubo) 1678 ihu=idbo(ncubo,i4) 1679 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 1680 end do 1681 end do 1682 endif 1683 1684 endif 1685 1686 end if 1687 end do 1688 1689 end if 1690 end do 1691 1692 40 continue 1693 end if 1694 1695 return 1696 end 1697********************************************************************** 1698********************************************************************** 1699 1700 subroutine molen 1701 1702********************************************************************** 1703#include "cbka.blk" 1704#include "cbkbo.blk" 1705#include "cbkconst.blk" 1706#include "cbkc.blk" 1707#include "cbkd.blk" 1708#include "cbkenergies.blk" 1709#include "cbkff.blk" 1710#include "cbkia.blk" 1711#include "cbkidbo.blk" 1712#include "cbknmolat.blk" 1713#include "cbknubon2.blk" 1714#include "control.blk" 1715#include "small.blk" 1716 dimension virial_tmp(3,3),virialsym(6) 1717********************************************************************** 1718* * 1719* Calculate molecular energy and first derivatives * 1720* Only used to prevent creating virtual electrons * 1721* * 1722********************************************************************** 1723c$$$ if (ndebug.eq.1) then 1724c$$$C open (65,file='fort.65',status='unknown',access='append') 1725c$$$ write (65,*) 'In molen' 1726c$$$ call timer(65) 1727c$$$ close (65) 1728c$$$ end if 1729 emol=zero 1730 return 1731 do i1=1,nmolo 1732 1733 enelm=0.0 1734 do i2=1,na 1735 if (ia(i2,3+mbond).eq.i1) then 1736 it1=ia(i2,1) 1737 enelm=enelm+aval(it1) 1738 end if 1739 end do 1740 1741 na1m=nmolat(i1,1) 1742 1743 enelm=2*int(enelm*0.50) 1744* enelm=elmol(i1) 1745 bomsum=zero 1746 do i2=1,na1m 1747 ihu=nmolat(i1,i2+1) 1748 do i3=1,ia(ihu,2) 1749 ihu2=nubon2(ihu,i3) 1750 bomsum=bomsum+bo(ihu2) 1751 end do 1752 end do 1753 diff=(bomsum-enelm) 1754 exphu=exp(-vpar(37)*diff) 1755 exphu2=1.0/(1.0+15.0*exphu) 1756 emolh=zero 1757 demoldsbo=zero 1758 emolh=vpar(38)*exphu2 1759 emol=emol+emolh 1760 demoldsbo=vpar(38)*vpar(37)*15.0*exphu2*exphu2*exphu 1761 1762 do i2=1,na1m 1763 ihu1=nmolat(i1,i2+1) 1764 do i3=1,ia(ihu1,2) 1765 iob=ia(ihu1,2+i3) 1766 ncubo=nubon2(ihu1,i3) 1767 1768 if (Lvirial.eq.1) then 1769 do k1=1,3 1770 do k2=1,3 1771 virial_tmp(k1,k2) = 0.0 1772 end do 1773 end do 1774 endif 1775 1776 do i4=1,idbo1(ncubo) 1777 ihu=idbo(ncubo,i4) 1778 do k1=1,3 1779 ftmp = demoldsbo*dbondc(ncubo,k1,i4) 1780 d(k1,ihu)=d(k1,ihu)+ftmp 1781 1782 if (Lvirial.eq.1) then 1783 do k1p=1,3 1784 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 1785 end do 1786 endif 1787 1788 end do 1789 end do 1790 if (Lvirial.eq.1) then 1791 virialsym(1) = virial_tmp(1,1) 1792 virialsym(2) = virial_tmp(2,2) 1793 virialsym(3) = virial_tmp(3,3) 1794 virialsym(4) = virial_tmp(1,2) 1795 virialsym(5) = virial_tmp(1,3) 1796 virialsym(6) = virial_tmp(2,3) 1797 do k1 = 1,6 1798 virial(k1) = virial(k1) + virialsym(k1) 1799 end do 1800 1801 if (Latomvirial.eq.1) then 1802 frac = 1.0d0/idbo1(ncubo) 1803 do k1 = 1,6 1804 vtmp = virialsym(k1)*frac 1805 do i4=1,idbo1(ncubo) 1806 ihu=idbo(ncubo,i4) 1807 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 1808 end do 1809 end do 1810 endif 1811 1812 endif 1813 1814 end do 1815 end do 1816 1817 1818 end do 1819 1820 1821 return 1822 end 1823********************************************************************** 1824********************************************************************** 1825 1826 subroutine valang 1827 1828********************************************************************** 1829#include "cbka.blk" 1830#include "cbkabo.blk" 1831#include "cbkbo.blk" 1832#include "cbkbopi.blk" 1833#include "cbkbopi2.blk" 1834#include "cbkconst.blk" 1835#include "cbkc.blk" 1836#include "cbkd.blk" 1837#include "cbkdbopi2ndc.blk" 1838#include "cbkdbopindc.blk" 1839#include "cbkdcell.blk" 1840#include "cbkdhdc.blk" 1841#include "cbkenergies.blk" 1842#include "cbkff.blk" 1843#include "cbkh.blk" 1844#include "cbkia.blk" 1845#include "cbkidbo.blk" 1846#include "cbklonpar.blk" 1847#include "cbknubon2.blk" 1848#include "cbkvalence.blk" 1849#include "control.blk" 1850#include "valang.blk" 1851#include "small.blk" 1852 dimension j(3) 1853 dimension virial_tmp(3,3),virialsym(6) 1854********************************************************************** 1855* * 1856* Calculate valency angle energies and first derivatives * 1857* * 1858********************************************************************** 1859c$$$ if (ndebug.eq.1) then 1860c$$$C open (65,file='fort.65',status='unknown',access='append') 1861c$$$ write (65,*) 'In valang' 1862c$$$ call timer(65) 1863c$$$ close (65) 1864c$$$ end if 1865* eco=0.0 1866 ev=0.0 1867 ecoa=0.0 1868 if (nval.eq.0) return 1869 1870 do 10 i1=1,nval 1871 ity=iv(i1,1) 1872 j(1)=iv(i1,2) 1873 j(2)=iv(i1,3) 1874 j(3)=iv(i1,4) 1875 1876 if (j(2) .le. na_local) then 1877 1878 la=iv(i1,5) 1879 lb=iv(i1,6) 1880 boa=bo(la)-cutof2 1881 bob=bo(lb)-cutof2 1882 if (boa.lt.zero.or.bob.lt.zero) goto 10 1883 1884 hl=h(i1) ! Calculated earlier in routine calval 1885********************************************************************** 1886* * 1887* Calculate valency angle energy * 1888* * 1889********************************************************************** 1890 nbocen=ia(j(2),2) 1891 sbo2=0.0 1892 vmbo=1.0 1893 1894 do i2=1,nbocen 1895 ibv=nubon2(j(2),i2) 1896 if (bo(ibv).gt.0.0) then 1897 vmbo=vmbo*exp(-bo(ibv)**8) 1898 sbo2=sbo2+bopi(ibv)+bopi2(ibv) 1899 endif 1900 end do 1901 1902 ity2=ia(j(2),1) 1903* exbo=abo(j(2))-stlp(ia(j(2),1)) 1904 exbo=abo(j(2))-valf(ity2) 1905* if (exbo.gt.zero) exbo=zero 1906* expov=exp(vka8(ity)*exbo) 1907* expov2=exp(-vpar(13)*exbo) 1908* htov1=2.0+expov2 1909* htov2=1.0+expov+expov2 1910* evboadj=htov1/htov2 1911 evboadj=1.0 1912 expun=exp(-vkac(ity)*exbo) 1913 expun2=exp(vpar(15)*exbo) 1914 htun1=2.0+expun2 1915 htun2=1.0+expun+expun2 1916 evboadj2=vval4(ity2)-(vval4(ity2)-1.0)*htun1/htun2 1917********************************************************************** 1918* * 1919* Calculate number of lone pairs * 1920* * 1921********************************************************************** 1922 dsbo2dvlp=(1.0-vmbo) 1923 vlpadj=zero 1924 exlp1=abo(j(2))-stlp(ia(j(2),1)) 1925 exlp2=2.0*int(exlp1/2.0) 1926 exlp=exlp1-exlp2 1927 if (exlp.lt.zero) then 1928* expvlp=exp(-vpar(16)*(2.0+exlp)*(2.0+exlp)) 1929* vlpadj=expvlp-int(exlp1/2.0) 1930* dsbo2dvlp=(1.0-vmbo)*(1.0-vpar(34)*2.0* 1931* $(2.0+exlp)*vpar(16)*expvlp) 1932 vlpadj=vlp(j(2)) 1933 dsbo2dvlp=(1.0-vmbo)*(1.0+vpar(34)*dvlpdsbo(j(2))) 1934 end if 1935 1936 sbo2=sbo2+(1.0-vmbo)*(-exbo-vpar(34)*vlpadj) 1937 dsbo2dvmbo=exbo+vpar(34)*vlpadj 1938 1939 sbo2h=sbo2 1940 powv=vpar(17) 1941 if (sbo2.le.0.0) sbo2h=0.0 1942 if (sbo2.gt.0.0.and.sbo2.le.1.0) sbo2h=sbo2**powv 1943 if (sbo2.gt.1.0.and.sbo2.lt.2.0) sbo2h=2.0-(2.0-sbo2)**powv 1944 if (sbo2.gt.2.0) sbo2h=2.0 1945 thba=th0(ity) 1946 expsbo=exp(-vpar(18)*(2.0-sbo2h)) 1947 thetao=180.0-thba*(1.0-expsbo) 1948 1949 thetao=thetao*dgrrdn 1950 thdif=(thetao-hl) 1951 thdi2=thdif*thdif 1952 dthsbo=dgrrdn*thba*vpar(18)*expsbo 1953 if (sbo2.lt.0.0) dthsbo=zero 1954 if (sbo2.gt.0.0.and.sbo2.le.1.0) 1955 $dthsbo=powv*(sbo2**(powv-1.0))*dgrrdn*thba*vpar(18)*expsbo 1956 if (sbo2.gt.1.0.and.sbo2.lt.2.0) 1957 $dthsbo=powv*((2.0-sbo2)**(powv-1.0))*dgrrdn*thba*vpar(18)*expsbo 1958 if (sbo2.gt.2.0) dthsbo=zero 1959 1960 exphu=vka(ity)*exp(-vka3(ity)*thdi2) 1961 exphu2=vka(ity)-exphu 1962 if (vka(ity).lt.zero) exphu2=exphu2-vka(ity) !To avoid linear Me-H-Me angles (6/6/06) 1963 boap=boa**vval2(ity) 1964 boap2=boa**(vval2(ity)-1.0) 1965 bobp=bob**vval2(ity) 1966 bobp2=bob**(vval2(ity)-1.0) 1967 exa=exp(-vval1(ity2)*boap) 1968 exb=exp(-vval1(ity2)*bobp) 1969 dexadboa=vval2(ity)*vval1(ity2)*exa*boap2 1970 dexbdbob=vval2(ity)*vval1(ity2)*exb*bobp2 1971 exa2=(1.0-exa) 1972 exb2=(1.0-exb) 1973 1974 evh=evboadj2*evboadj*exa2*exb2*exphu2 1975 devdlb=evboadj2*evboadj*dexbdbob*exa2*exphu2 1976 devdla=evboadj2*evboadj*dexadboa*exb2*exphu2 1977 devdsbo=2.0*evboadj2*evboadj*dthsbo*exa2*exb2* 1978 $vka3(ity)*thdif*exphu 1979 devdh=-2.0*evboadj2*evboadj*exa2*exb2*vka3(ity)*thdif*exphu 1980 1981 devdsbo2= 1982 $evboadj*exa2*exb2*exphu2*(vval4(ity2)-1.0)*(-vpar(15)*expun2/htun2 1983 $+htun1*(vpar(15)*expun2-vkac(ity)*expun)/(htun2*htun2)) 1984 1985* devdsbo2=-evboadj2*exa2*exb2*exphu2*(vpar(13)*expov2/htov2+ 1986* $htov1*(vka8(ity)*expov-vpar(13)*expov2)/(htov2*htov2))+ 1987* $evboadj*exa2*exb2*exphu2*(vpar(14)-1.0)*(-vpar(15)*expun2/htun2 1988* $+htun1*(vpar(15)*expun2-vkac(ity)*expun)/(htun2*htun2)) 1989 1990 if (j(2) .le. na_local) then 1991 ev=ev+evh 1992 estrain(j(2))=estrain(j(2))+evh !central atom energy 1993 endif 1994 1995* write (64,'(4i8,18f8.2)')mdstep,j(1),j(2),j(3),sbo2,sbo2h, 1996* $thetao*rdndgr,hl*rdndgr,bo(la),bo(lb),bopi(la), 1997* $vlp(j(2)),exbo,vlpadj,vmbo,evh,ev,vka(ity) 1998********************************************************************** 1999* * 2000* Calculate penalty for two double bonds in valency angle * 2001* * 2002********************************************************************** 2003 exbo=abo(j(2))-aval(ia(j(2),1)) 2004 expov=exp(vpar(22)*exbo) 2005 expov2=exp(-vpar(21)*exbo) 2006 htov1=2.0+expov2 2007 htov2=1.0+expov+expov2 2008 ecsboadj=htov1/htov2 2009 exphu1=exp(-vpar(20)*(boa-2.0)*(boa-2.0)) 2010 exphu2=exp(-vpar(20)*(bob-2.0)*(bob-2.0)) 2011 2012 epenh=vkap(ity)*ecsboadj*exphu1*exphu2 2013 estrain(j(2))=estrain(j(2))+epenh 2014 epen=epen+epenh 2015 decoadboa=-2.0*vpar(20)*epenh*(boa-2.0) 2016 decoadbob=-2.0*vpar(20)*epenh*(bob-2.0) 2017 2018 decdsbo2=-vkap(ity)*exphu1*exphu2*(vpar(21)*expov2/htov2+htov1* 2019 $(vpar(22)*expov-vpar(21)*expov2)/(htov2*htov2)) 2020********************************************************************** 2021* * 2022* Calculate valency angle conjugation energy * 2023* * 2024********************************************************************** 2025 unda=abo(j(1))-boa 2026* ovb=abo(j(2))-valf(ia(j(2),1)) 2027 ovb=abo(j(2))-vval3(ia(j(2),1)) !Modification for Ru 7/6/2004 2028 2029 undc=abo(j(3))-bob 2030 ba=(boa-1.50)*(boa-1.50) 2031 bb=(bob-1.50)*(bob-1.50) 2032 exphua=exp(-vpar(31)*ba) 2033 exphub=exp(-vpar(31)*bb) 2034 exphuua=exp(-vpar(39)*unda*unda) 2035 exphuob=exp(vpar(3)*ovb) 2036 exphuuc=exp(-vpar(39)*undc*undc) 2037 hulpob=1.0/(1.0+exphuob) 2038 ecoah=vka8(ity)*exphua*exphub*exphuua*exphuuc*hulpob 2039 decodbola=-2.0*vka8(ity)*(boa-1.50)*vpar(31)*exphua*exphub 2040 $*exphuua*exphuuc*hulpob+vpar(39)*vka8(ity)*exphua*exphub* 2041 $exphuua*exphuuc*hulpob*2.0*unda 2042 decodbolb=-2.0*vka8(ity)*(bob-1.50)*vpar(31)*exphua*exphub 2043 $*exphuua*exphuuc*hulpob+vpar(39)*vka8(ity)*exphua*exphub* 2044 $exphuua*exphuuc*hulpob*2.0*undc 2045 decodboua=-2.0*unda*vka8(ity)*vpar(39)*exphua*exphub 2046 $*exphuua*exphuuc*hulpob 2047 decodbouc=-2.0*undc*vka8(ity)*vpar(39)*exphua*exphub 2048 $*exphuua*exphuuc*hulpob 2049 decodboob=-vka8(ity)*exphua*exphub*exphuua*exphuuc*hulpob* 2050 $hulpob*vpar(3)*exphuob 2051* decodboob=zero 2052* decodboua=zero 2053* decodbouc=zero 2054 2055 ecoa=ecoa+ecoah 2056 estrain(j(2))=estrain(j(2))+ecoah !central atom energy 2057 2058********************************************************************** 2059* * 2060* Calculate derivative valency energy to cartesian coordinates * 2061* * 2062********************************************************************** 2063 if (Lvirial.eq.1) then 2064 do k1=1,3 2065 do k2=1,3 2066 virial_tmp(k1,k2) = 0.0 2067 end do 2068 end do 2069 endif 2070 2071 do k1=1,3 2072 do k2=1,3 2073 ftmp = devdh*dhdc(i1,k1,k2) 2074 d(k1,j(k2))=d(k1,j(k2))+ftmp 2075 2076 if (Lvirial.eq.1) then 2077 do k1p=1,3 2078 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(j(k2),k1p) 2079 end do 2080 endif 2081 2082 end do 2083 end do 2084 if (Lvirial.eq.1) then 2085 virialsym(1) = virial_tmp(1,1) 2086 virialsym(2) = virial_tmp(2,2) 2087 virialsym(3) = virial_tmp(3,3) 2088 virialsym(4) = virial_tmp(1,2) 2089 virialsym(5) = virial_tmp(1,3) 2090 virialsym(6) = virial_tmp(2,3) 2091 do k1 = 1,6 2092 virial(k1) = virial(k1) + virialsym(k1) 2093 end do 2094 2095 if (Latomvirial.eq.1) then 2096 frac = 1.0d0/3 2097 do k1 = 1,6 2098 vtmp = virialsym(k1)*frac 2099 do k2=1,3 2100 ihu=j(k2) 2101 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 2102 end do 2103 end do 2104 endif 2105 2106 endif 2107 2108 if (Lvirial.eq.1) then 2109 do k1=1,3 2110 do k2=1,3 2111 virial_tmp(k1,k2) = 0.0 2112 end do 2113 end do 2114 endif 2115 2116 do i2=1,idbo1(la) 2117 ihu=idbo(la,i2) 2118 do k1=1,3 2119 ftmp = (devdla+decoadboa+decodbola)* 2120 $dbondc(la,k1,i2) 2121 d(k1,ihu)=d(k1,ihu)+ftmp 2122 2123 if (Lvirial.eq.1) then 2124 do k1p=1,3 2125 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 2126 end do 2127 endif 2128 2129 end do 2130 end do 2131 if (Lvirial.eq.1) then 2132 virialsym(1) = virial_tmp(1,1) 2133 virialsym(2) = virial_tmp(2,2) 2134 virialsym(3) = virial_tmp(3,3) 2135 virialsym(4) = virial_tmp(1,2) 2136 virialsym(5) = virial_tmp(1,3) 2137 virialsym(6) = virial_tmp(2,3) 2138 do k1 = 1,6 2139 virial(k1) = virial(k1) + virialsym(k1) 2140 end do 2141 2142 if (Latomvirial.eq.1) then 2143 frac = 1.0d0/idbo1(la) 2144 do k1 = 1,6 2145 vtmp = virialsym(k1)*frac 2146 do i2=1,idbo1(la) 2147 ihu=idbo(la,i2) 2148 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 2149 end do 2150 end do 2151 endif 2152 2153 endif 2154 2155 if (Lvirial.eq.1) then 2156 do k1=1,3 2157 do k2=1,3 2158 virial_tmp(k1,k2) = 0.0 2159 end do 2160 end do 2161 endif 2162 2163 do i2=1,idbo1(lb) 2164 ihu=idbo(lb,i2) 2165 do k1=1,3 2166 ftmp = (devdlb+decoadbob+decodbolb)* 2167 $dbondc(lb,k1,i2) 2168 d(k1,ihu)=d(k1,ihu)+ftmp 2169 2170 if (Lvirial.eq.1) then 2171 do k1p=1,3 2172 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 2173 end do 2174 endif 2175 2176 end do 2177 end do 2178 if (Lvirial.eq.1) then 2179 virialsym(1) = virial_tmp(1,1) 2180 virialsym(2) = virial_tmp(2,2) 2181 virialsym(3) = virial_tmp(3,3) 2182 virialsym(4) = virial_tmp(1,2) 2183 virialsym(5) = virial_tmp(1,3) 2184 virialsym(6) = virial_tmp(2,3) 2185 do k1 = 1,6 2186 virial(k1) = virial(k1) + virialsym(k1) 2187 end do 2188 2189 if (Latomvirial.eq.1) then 2190 frac = 1.0d0/idbo1(lb) 2191 do k1 = 1,6 2192 vtmp = virialsym(k1)*frac 2193 do i2=1,idbo1(lb) 2194 ihu=idbo(lb,i2) 2195 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 2196 end do 2197 end do 2198 endif 2199 2200 endif 2201 2202 do i2=1,nbocen 2203 j5=ia(j(2),2+i2) 2204 ibv=nubon2(j(2),i2) 2205 if (bo(ibv).gt.0.0) then 2206 dvmbodbo=-vmbo*8.0*bo(ibv)**7 2207 2208 if (Lvirial.eq.1) then 2209 do k1=1,3 2210 do k2=1,3 2211 virial_tmp(k1,k2) = 0.0 2212 end do 2213 end do 2214 endif 2215 2216 do i3=1,idbo1(ibv) 2217 ihu=idbo(ibv,i3) 2218 do k1=1,3 2219 ftmp = (-dsbo2dvlp*devdsbo+devdsbo2+decdsbo2 2220 $+dvmbodbo*dsbo2dvmbo*devdsbo)* 2221 $dbondc(ibv,k1,i3)+devdsbo*(dbopindc(ibv,k1,i3)+ 2222 $dbopi2ndc(ibv,k1,i3)) 2223 d(k1,ihu)=d(k1,ihu)+ftmp 2224 2225 if (Lvirial.eq.1) then 2226 do k1p=1,3 2227 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 2228 end do 2229 endif 2230 2231 end do 2232 end do 2233 if (Lvirial.eq.1) then 2234 virialsym(1) = virial_tmp(1,1) 2235 virialsym(2) = virial_tmp(2,2) 2236 virialsym(3) = virial_tmp(3,3) 2237 virialsym(4) = virial_tmp(1,2) 2238 virialsym(5) = virial_tmp(1,3) 2239 virialsym(6) = virial_tmp(2,3) 2240 do k1 = 1,6 2241 virial(k1) = virial(k1) + virialsym(k1) 2242 end do 2243 2244 if (Latomvirial.eq.1) then 2245 frac = 1.0d0/idbo1(ibv) 2246 do k1 = 1,6 2247 vtmp = virialsym(k1)*frac 2248 do i3=1,idbo1(ibv) 2249 ihu=idbo(ibv,i3) 2250 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 2251 end do 2252 end do 2253 endif 2254 2255 endif 2256 2257 endif 2258 end do 2259 2260 do i2=1,ia(j(1),2) 2261 j5=ia(j(1),2+i2) 2262 ibv=nubon2(j(1),i2) 2263 if (bo(ibv).gt.0.0) then 2264 2265 if (Lvirial.eq.1) then 2266 do k1=1,3 2267 do k2=1,3 2268 virial_tmp(k1,k2) = 0.0 2269 end do 2270 end do 2271 endif 2272 2273 do i3=1,idbo1(ibv) 2274 ihu=idbo(ibv,i3) 2275 do k1=1,3 2276 ftmp = decodboua*dbondc(ibv,k1,i3) 2277 d(k1,ihu)=d(k1,ihu)+ftmp 2278 2279 if (Lvirial.eq.1) then 2280 do k1p=1,3 2281 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 2282 end do 2283 endif 2284 2285 end do 2286 end do 2287 if (Lvirial.eq.1) then 2288 virialsym(1) = virial_tmp(1,1) 2289 virialsym(2) = virial_tmp(2,2) 2290 virialsym(3) = virial_tmp(3,3) 2291 virialsym(4) = virial_tmp(1,2) 2292 virialsym(5) = virial_tmp(1,3) 2293 virialsym(6) = virial_tmp(2,3) 2294 do k1 = 1,6 2295 virial(k1) = virial(k1) + virialsym(k1) 2296 end do 2297 2298 if (Latomvirial.eq.1) then 2299 frac = 1.0d0/idbo1(ibv) 2300 do k1 = 1,6 2301 vtmp = virialsym(k1)*frac 2302 do i3=1,idbo1(ibv) 2303 ihu=idbo(ibv,i3) 2304 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 2305 end do 2306 end do 2307 endif 2308 2309 endif 2310 2311 endif 2312 end do 2313 2314 do i2=1,ia(j(2),2) 2315 j5=ia(j(2),2+i2) 2316 ibv=nubon2(j(2),i2) 2317 if (bo(ibv).gt.0.0) then 2318 2319 if (Lvirial.eq.1) then 2320 do k1=1,3 2321 do k2=1,3 2322 virial_tmp(k1,k2) = 0.0 2323 end do 2324 end do 2325 endif 2326 2327 do i3=1,idbo1(ibv) 2328 ihu=idbo(ibv,i3) 2329 do k1=1,3 2330 ftmp = decodboob*dbondc(ibv,k1,i3) 2331 d(k1,ihu)=d(k1,ihu)+ftmp 2332 2333 if (Lvirial.eq.1) then 2334 do k1p=1,3 2335 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 2336 end do 2337 endif 2338 2339 end do 2340 end do 2341 if (Lvirial.eq.1) then 2342 virialsym(1) = virial_tmp(1,1) 2343 virialsym(2) = virial_tmp(2,2) 2344 virialsym(3) = virial_tmp(3,3) 2345 virialsym(4) = virial_tmp(1,2) 2346 virialsym(5) = virial_tmp(1,3) 2347 virialsym(6) = virial_tmp(2,3) 2348 do k1 = 1,6 2349 virial(k1) = virial(k1) + virialsym(k1) 2350 end do 2351 2352 if (Latomvirial.eq.1) then 2353 frac = 1.0d0/idbo1(ibv) 2354 do k1 = 1,6 2355 vtmp = virialsym(k1)*frac 2356 do i3=1,idbo1(ibv) 2357 ihu=idbo(ibv,i3) 2358 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 2359 end do 2360 end do 2361 endif 2362 2363 endif 2364 2365 endif 2366 end do 2367 2368 do i2=1,ia(j(3),2) 2369 j5=ia(j(3),2+i2) 2370 ibv=nubon2(j(3),i2) 2371 if (bo(ibv).gt.0.0) then 2372 2373 if (Lvirial.eq.1) then 2374 do k1=1,3 2375 do k2=1,3 2376 virial_tmp(k1,k2) = 0.0 2377 end do 2378 end do 2379 endif 2380 2381 do i3=1,idbo1(ibv) 2382 ihu=idbo(ibv,i3) 2383 do k1=1,3 2384 ftmp = decodbouc*dbondc(ibv,k1,i3) 2385 d(k1,ihu)=d(k1,ihu)+ftmp 2386 2387 if (Lvirial.eq.1) then 2388 do k1p=1,3 2389 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 2390 end do 2391 endif 2392 2393 end do 2394 end do 2395 if (Lvirial.eq.1) then 2396 virialsym(1) = virial_tmp(1,1) 2397 virialsym(2) = virial_tmp(2,2) 2398 virialsym(3) = virial_tmp(3,3) 2399 virialsym(4) = virial_tmp(1,2) 2400 virialsym(5) = virial_tmp(1,3) 2401 virialsym(6) = virial_tmp(2,3) 2402 do k1 = 1,6 2403 virial(k1) = virial(k1) + virialsym(k1) 2404 end do 2405 2406 if (Latomvirial.eq.1) then 2407 frac = 1.0d0/idbo1(ibv) 2408 do k1 = 1,6 2409 vtmp = virialsym(k1)*frac 2410 do i3=1,idbo1(ibv) 2411 ihu=idbo(ibv,i3) 2412 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 2413 end do 2414 end do 2415 endif 2416 2417 endif 2418 2419 endif 2420 end do 2421 2422 endif 2423 2424 10 continue 2425 2426 return 2427 end 2428********************************************************************** 2429********************************************************************** 2430 2431 subroutine hbond 2432 2433********************************************************************** 2434#include "cbka.blk" 2435#include "cbkbo.blk" 2436#include "cbkconst.blk" 2437#include "cbkc.blk" 2438#include "cbkd.blk" 2439#include "cbkdcell.blk" 2440#include "cbkenergies.blk" 2441#include "cbkidbo.blk" 2442#include "cbksrthb.blk" 2443#include "control.blk" 2444#include "cbkhbond.blk" 2445#include "small.blk" 2446 dimension drda(3),j(3),dvdc(3,3),dargdc(3,3) 2447 dimension virial_tmp(3,3),virialsym(6) 2448********************************************************************** 2449* * 2450* Calculate hydrogen bond energies and first derivatives * 2451* * 2452********************************************************************** 2453c$$$ if (ndebug.eq.1) then 2454c$$$C open (65,file='fort.65',status='unknown',access='append') 2455c$$$ write (65,*) 'In hbond' 2456c$$$ call timer(65) 2457c$$$ close (65) 2458c$$$ end if 2459 ehb=zero 2460 do 10 i1=1,nhb 2461 ityhb=ihb(i1,1) 2462 j(1)=ihb(i1,2) 2463 j(2)=ihb(i1,3) 2464 j(3)=ihb(i1,4) 2465 la=ihb(i1,5) 2466 boa=bo(la) 2467 call dista2(j(2),j(3),rda,dxm,dym,dzm) 2468 drda(1)=dxm/rda 2469 drda(2)=dym/rda 2470 drda(3)=dzm/rda 2471 call calvalhb(j(1),j(2),j(3),ix,iy,iz,arg,hhb(i1),dvdc,dargdc) 2472 rhu1=rhb(ityhb)/rda 2473 rhu2=rda/rhb(ityhb) 2474 sinhu=sin(hhb(i1)/2.0) 2475 sin2=sinhu*sinhu 2476 exphu1=exp(-vhb1(ityhb)*boa) 2477 exphu2=exp(-vhb2(ityhb)*(rhu1+rhu2-2.0)) 2478 if (lhbnew .eq. 0) then 2479 ehbh=(1.0-exphu1)*dehb(ityhb)*exphu2*sin2*sin2*sin2*sin2 2480 else 2481 ehbh=(1.0-exphu1)*dehb(ityhb)*exphu2*sin2*sin2 2482 endif 2483 ehb=ehb+ehbh 2484 estrain(j(2))=estrain(j(2))+ehbh !2nd atom energy 2485 2486********************************************************************** 2487* * 2488* Calculate first derivatives * 2489* * 2490********************************************************************** 2491 if (lhbnew .eq. 0) then 2492 dehbdbo=vhb1(ityhb)*exphu1*dehb(ityhb)*exphu2*sin2*sin2* 2493 $ sin2*sin2 2494 dehbdv=(1.0-exphu1)*dehb(ityhb)*exphu2* 2495 $ 4.0*sin2*sin2*sin2*sinhu*cos(hhb(i1)/2.0) 2496 dehbdrda=(1.0-exphu1)*dehb(ityhb)*sin2*sin2*sin2*sin2* 2497 $ vhb2(ityhb)*(rhb(ityhb)/(rda*rda)-1.0/rhb(ityhb))*exphu2 2498 else 2499 dehbdbo=vhb1(ityhb)*exphu1*dehb(ityhb)*exphu2*sin2*sin2 2500 dehbdv=(1.0-exphu1)*dehb(ityhb)*exphu2* 2501 $ 2.0*sin2*sinhu*cos(hhb(i1)/2.0) 2502 dehbdrda=(1.0-exphu1)*dehb(ityhb)*sin2*sin2* 2503 $ vhb2(ityhb)*(rhb(ityhb)/(rda*rda)-1.0/rhb(ityhb))*exphu2 2504 endif 2505 2506 if (Lvirial.eq.1) then 2507 do k1=1,3 2508 do k2=1,3 2509 virial_tmp(k1,k2) = 0.0 2510 end do 2511 end do 2512 endif 2513 2514 do k1=1,3 2515 ftmp = dehbdrda*drda(k1) 2516 d(k1,j(2))=d(k1,j(2))+ftmp 2517 d(k1,j(3))=d(k1,j(3))-ftmp 2518 2519 if (Lvirial.eq.1) then 2520 do k1p=1,3 2521 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ 2522 $ ftmp*c(j(2),k1p)-ftmp*c(j(3),k1p) 2523 end do 2524 endif 2525 2526 end do 2527 if (Lvirial.eq.1) then 2528 virialsym(1) = virial_tmp(1,1) 2529 virialsym(2) = virial_tmp(2,2) 2530 virialsym(3) = virial_tmp(3,3) 2531 virialsym(4) = virial_tmp(1,2) 2532 virialsym(5) = virial_tmp(1,3) 2533 virialsym(6) = virial_tmp(2,3) 2534 do k1 = 1,6 2535 virial(k1) = virial(k1) + virialsym(k1) 2536 end do 2537 2538 if (Latomvirial.eq.1) then 2539 frac = 1.0d0/2 2540 do k1 = 1,6 2541 vtmp = virialsym(k1)*frac 2542 ihu = j(2) 2543 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 2544 ihu = j(3) 2545 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 2546 end do 2547 endif 2548 2549 endif 2550 2551 if (Lvirial.eq.1) then 2552 do k1=1,3 2553 do k2=1,3 2554 virial_tmp(k1,k2) = 0.0 2555 end do 2556 end do 2557 endif 2558 2559 do k1=1,3 2560 do k2=1,3 2561 ftmp = dehbdv*dvdc(k1,k2) 2562 d(k1,j(k2))=d(k1,j(k2))+ftmp 2563 2564 if (Lvirial.eq.1) then 2565 do k1p=1,3 2566 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(j(k2),k1p) 2567 end do 2568 endif 2569 2570 end do 2571 end do 2572 if (Lvirial.eq.1) then 2573 virialsym(1) = virial_tmp(1,1) 2574 virialsym(2) = virial_tmp(2,2) 2575 virialsym(3) = virial_tmp(3,3) 2576 virialsym(4) = virial_tmp(1,2) 2577 virialsym(5) = virial_tmp(1,3) 2578 virialsym(6) = virial_tmp(2,3) 2579 do k1 = 1,6 2580 virial(k1) = virial(k1) + virialsym(k1) 2581 end do 2582 2583 if (Latomvirial.eq.1) then 2584 frac = 1.0d0/3 2585 do k1 = 1,6 2586 vtmp = virialsym(k1)*frac 2587 do k2=1,3 2588 ihu=j(k2) 2589 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 2590 end do 2591 end do 2592 endif 2593 2594 endif 2595 2596 if (Lvirial.eq.1) then 2597 do k1=1,3 2598 do k2=1,3 2599 virial_tmp(k1,k2) = 0.0 2600 end do 2601 end do 2602 endif 2603 2604 do i2=1,idbo1(la) 2605 ihu=idbo(la,i2) 2606 do k1=1,3 2607 ftmp = dehbdbo*dbondc(la,k1,i2) 2608 d(k1,ihu)=d(k1,ihu)+ftmp 2609 2610 if (Lvirial.eq.1) then 2611 do k1p=1,3 2612 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 2613 end do 2614 endif 2615 2616 end do 2617 end do 2618 if (Lvirial.eq.1) then 2619 virialsym(1) = virial_tmp(1,1) 2620 virialsym(2) = virial_tmp(2,2) 2621 virialsym(3) = virial_tmp(3,3) 2622 virialsym(4) = virial_tmp(1,2) 2623 virialsym(5) = virial_tmp(1,3) 2624 virialsym(6) = virial_tmp(2,3) 2625 do k1 = 1,6 2626 virial(k1) = virial(k1) + virialsym(k1) 2627 end do 2628 2629 if (Latomvirial.eq.1) then 2630 frac = 1.0d0/idbo1(la) 2631 do k1 = 1,6 2632 vtmp = virialsym(k1)*frac 2633 do i2=1,idbo1(la) 2634 ihu=idbo(la,i2) 2635 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 2636 end do 2637 end do 2638 endif 2639 2640 endif 2641 2642 10 continue 2643 return 2644 end 2645 2646********************************************************************** 2647********************************************************************** 2648 2649 subroutine torang 2650 2651********************************************************************** 2652#include "cbka.blk" 2653#include "cbkabo.blk" 2654#include "cbkbo.blk" 2655#include "cbkbopi.blk" 2656#include "cbkc.blk" 2657#include "cbkconst.blk" 2658#include "cbkd.blk" 2659#include "cbkdbopindc.blk" 2660#include "cbkdcell.blk" 2661#include "cbkdhdc.blk" 2662#include "cbkdrdc.blk" 2663#include "cbkenergies.blk" 2664#include "cbkff.blk" 2665#include "cbkfftorang.blk" 2666#include "cbkh.blk" 2667#include "cbkia.blk" 2668#include "cbkidbo.blk" 2669#include "cbkinit.blk" 2670#include "cbknubon2.blk" 2671#include "cbkrbo.blk" 2672#include "cbktorang.blk" 2673#include "cbktorsion.blk" 2674#include "cbktregime.blk" 2675#include "cbkvalence.blk" 2676#include "cellcoord.blk" 2677#include "control.blk" 2678#include "small.blk" 2679 2680 DIMENSION A(3),DRDA(3),DADC(4),DRADC(3,4),DRBDC(3,4), 2681 $DRCDC(3,4),DHDDC(3,4),DHEDC(3,4),DRVDC(3,4),DTDC(3,4), 2682 $DNDC(3,4) 2683 dimension j(4),dh1rdc(3,3),dh2rdc(3,3),dargdc(3,3) 2684 dimension virial_tmp(3,3),virialsym(6) 2685********************************************************************** 2686* * 2687* Calculate torsion angle energies and first derivatives * 2688* * 2689********************************************************************** 2690c$$$ if (ndebug.eq.1) then 2691c$$$C open (65,file='fort.65',status='unknown',access='append') 2692c$$$ write (65,*) 'In torang' 2693c$$$ call timer(65) 2694c$$$ close (65) 2695c$$$ end if 2696 do k1=1,3 2697 do k2=1,4 2698 dhddc(k1,k2)=0.0 2699 dhedc(k1,k2)=0.0 2700 dradc(k1,k2)=0.0 2701 drbdc(k1,k2)=0.0 2702 drcdc(k1,k2)=0.0 2703 end do 2704 end do 2705 et=0.0 2706 eth12=0.0 2707 eco=0.0 2708 dadc(1)=1.0 2709 dadc(2)=0.0 2710 dadc(3)=0.0 2711 dadc(4)=-1.0 2712 if (ntor.eq.0) return 2713 2714 do 10 i1=1,ntor 2715 j(1)=it(i1,2) 2716 j(2)=it(i1,3) 2717 j(3)=it(i1,4) 2718 j(4)=it(i1,5) 2719 2720 ity=it(i1,1) 2721 la=it(i1,6) 2722 lb=it(i1,7) 2723 lc=it(i1,8) 2724 call calvalres(j(1),j(2),j(3),arg1,ht1,dh1rdc,dargdc) 2725 call calvalres(j(2),j(3),j(4),arg2,ht2,dh2rdc,dargdc) 2726 boa=bo(la)-cutof2 2727 bob=bo(lb)-cutof2 2728 boc=bo(lc)-cutof2 2729 if (boa.lt.zero.or.bob.lt.zero.or.boc.lt.zero) 2730 $goto 10 2731 r42=0.0 2732 ivl1=ibsym(la) 2733 ivl2=ibsym(lb) 2734 ivl3=ibsym(lc) 2735 isign1=1 2736 isign2=1 2737 isign3=1 2738 rla=rbo(la) 2739 rlb=rbo(lb) 2740 2741 call dista2(j(1),j(4),r4,a(1),a(2),a(3)) 2742********************************************************************** 2743* * 2744* Determine torsion angle * 2745* * 2746********************************************************************** 2747 d142=r4*r4 2748 rla=rbo(la) 2749 rlb=rbo(lb) 2750 rlc=rbo(lc) 2751 coshd=cos(ht1) 2752 coshe=cos(ht2) 2753 sinhd=sin(ht1) 2754 sinhe=sin(ht2) 2755 poem=2.0*rla*rlc*sinhd*sinhe 2756 poem2=poem*poem 2757 tel=rla*rla+rlb*rlb+rlc*rlc-d142-2.0*(rla*rlb*coshd-rla*rlc* 2758 $coshd*coshe+rlb*rlc*coshe) 2759 if (poem.lt.1e-20) poem=1e-20 2760 arg=tel/poem 2761 if (arg.gt.1.0) arg=1.0 2762 if (arg.lt.-1.0) arg=-1.0 2763 arg2=arg*arg 2764 thg(i1)=acos(arg)*rdndgr 2765 k1=j(1) 2766 k2=j(2) 2767 k3=j(3) 2768 k4=j(4) 2769 call dista2(k3,k2,dis,x3,y3,z3) 2770 y32z32=y3*y3+z3*z3 2771 wort1=sqrt(y32z32)+1e-6 2772 wort2=sqrt(y32z32+x3*x3)+1e-6 2773* if (wort1.lt.1e-6) wort1=1e-6 2774* if (wort2.lt.1e-6) wort2=1e-6 2775 sinalf=y3/wort1 2776 cosalf=z3/wort1 2777 sinbet=x3/wort2 2778 cosbet=wort1/wort2 2779 call dista2(k1,k2,dis,x1,y1,z1) 2780 x1=x1*cosbet-y1*sinalf*sinbet-z1*cosalf*sinbet 2781 y1=y1*cosalf-z1*sinalf 2782 wort3=sqrt(x1*x1+y1*y1)+1e-6 2783* if (wort3.lt.1e-6) wort3=1e-6 2784 singam=y1/wort3 2785 cosgam=x1/wort3 2786 call dista2(k4,k2,dis,x4,y4,z4) 2787 x4=x4*cosbet-y4*sinalf*sinbet-z4*cosalf*sinbet 2788 y4=y4*cosalf-z4*sinalf 2789 y4=x4*singam-y4*cosgam 2790 if (y4.gt.0.0) thg(i1)=-thg(i1) 2791 if (thg(i1).lt.-179.999999d0) thg(i1)=-179.999999d0 2792 if (thg(i1).gt.179.999999d0) thg(i1)=179.999999d0 2793 th2=thg(i1)*dgrrdn 2794********************************************************************** 2795* * 2796* Calculate torsion angle energy * 2797* * 2798********************************************************************** 2799 exbo1=abo(j(2))-valf(ia(j(2),1)) 2800 exbo2=abo(j(3))-valf(ia(j(3),1)) 2801 htovt=exbo1+exbo2 2802 expov=exp(vpar(26)*htovt) 2803 expov2=exp(-vpar(25)*(htovt)) 2804 htov1=2.0+expov2 2805 htov2=1.0+expov+expov2 2806 etboadj=htov1/htov2 2807 2808 btb2=bopi(lb)-1.0+etboadj 2809 bo2t=1.0-btb2 2810 bo2p=bo2t*bo2t 2811 bocor2=exp(v4(ity)*bo2p) 2812 2813 hsin=sinhd*sinhe 2814 ethhulp=0.50*v1(ity)*(1.0+arg)+v2(ity)*bocor2*(1.0-arg2)+ 2815 $v3(ity)*(0.50+2.0*arg2*arg-1.50*arg) 2816 2817 exphua=exp(-vpar(24)*boa) 2818 exphub=exp(-vpar(24)*bob) 2819 exphuc=exp(-vpar(24)*boc) 2820 bocor4=(1.0-exphua)*(1.0-exphub)*(1.0-exphuc) 2821 eth=hsin*ethhulp*bocor4 2822 2823 detdar=hsin*bocor4*(0.50*v1(ity)-2.0*v2(ity)*bocor2*arg+ 2824 $v3(ity)*(6.0*arg2-1.5d0)) 2825 detdhd=coshd*sinhe*bocor4*ethhulp 2826 detdhe=sinhd*coshe*bocor4*ethhulp 2827 2828 detdboa=vpar(24)*exphua*(1.0-exphub)*(1.0-exphuc)*ethhulp*hsin 2829 detdbopib=-bocor4*2.0*v4(ity)*v2(ity)* 2830 $bo2t*bocor2*(1.0-arg2)*hsin 2831 detdbob=vpar(24)*exphub*(1.0-exphua)* 2832 $(1.0-exphuc)*ethhulp*hsin 2833 detdboc=vpar(24)*exphuc*(1.0-exphua)* 2834 $(1.0-exphub)*ethhulp*hsin 2835 2836 detdsbo1=-(detdbopib)* 2837 $(vpar(25)*expov2/htov2+htov1* 2838 $(vpar(26)*expov-vpar(25)*expov2)/(htov2*htov2)) 2839 2840 et=et+eth 2841 estrain(j(2))=estrain(j(2))+0.50*eth !2nd atom energy 2842 estrain(j(3))=estrain(j(3))+0.50*eth !3rd atom energy 2843 2844********************************************************************** 2845* * 2846* Calculate conjugation energy * 2847* * 2848********************************************************************** 2849 ba=(boa-1.50)*(boa-1.50) 2850 bb=(bob-1.50)*(bob-1.50) 2851 bc=(boc-1.50)*(boc-1.50) 2852 exphua1=exp(-vpar(28)*ba) 2853 exphub1=exp(-vpar(28)*bb) 2854 exphuc1=exp(-vpar(28)*bc) 2855 sbo=exphua1*exphub1*exphuc1 2856 dbohua=-2.0*(boa-1.50)*vpar(28)*exphua1*exphub1*exphuc1 2857 dbohub=-2.0*(bob-1.50)*vpar(28)*exphua1*exphub1*exphuc1 2858 dbohuc=-2.0*(boc-1.50)*vpar(28)*exphua1*exphub1*exphuc1 2859 arghu0=(arg2-1.0)*sinhd*sinhe 2860 ehulp=vconj(ity)*(arghu0+1.0) 2861 ecoh=ehulp*sbo 2862 decodar=sbo*vconj(ity)*2.0*arg*sinhd*sinhe 2863 decodbola=dbohua*ehulp 2864 decodbolb=dbohub*ehulp 2865 decodbolc=dbohuc*ehulp 2866 decodhd=coshd*sinhe*vconj(ity)*sbo*(arg2-1.0) 2867 decodhe=coshe*sinhd*vconj(ity)*sbo*(arg2-1.0) 2868 2869 eco=eco+ecoh 2870 estrain(j(2))=estrain(j(2))+0.50*ecoh !2nd atom energy 2871 estrain(j(3))=estrain(j(3))+0.50*ecoh !3rd atom energy 2872 2873 1 continue 2874********************************************************************** 2875* * 2876* Calculate derivative torsion angle and conjugation energy * 2877* to cartesian coordinates * 2878* * 2879********************************************************************** 2880 SINTH=SIN(THG(i1)*DGRRDN) 2881 IF (SINTH.GE.0.0.AND.SINTH.LT.1.0D-10) SINTH=1.0D-10 2882 IF (SINTH.LT.0.0.AND.SINTH.GT.-1.0D-10) SINTH=-1.0D-10 2883 IF (j(1).EQ.IB(LA,2)) THEN 2884 DO K1=1,3 2885 DRADC(K1,1)=DRDC(LA,K1,1) 2886 DRADC(K1,2)=DRDC(LA,K1,2) 2887 end do 2888 ELSE 2889 DO K1=1,3 2890 DRADC(K1,1)=DRDC(LA,K1,2) 2891 DRADC(K1,2)=DRDC(LA,K1,1) 2892 end do 2893 ENDIF 2894 IF (j(2).EQ.IB(LB,2)) THEN 2895 DO K1=1,3 2896 DRBDC(K1,2)=DRDC(LB,K1,1) 2897 DRBDC(K1,3)=DRDC(LB,K1,2) 2898 end do 2899 ELSE 2900 DO K1=1,3 2901 DRBDC(K1,2)=DRDC(LB,K1,2) 2902 DRBDC(K1,3)=DRDC(LB,K1,1) 2903 end do 2904 ENDIF 2905 IF (j(3).EQ.IB(LC,2)) THEN 2906 DO K1=1,3 2907 DRCDC(K1,3)=DRDC(LC,K1,1) 2908 DRCDC(K1,4)=DRDC(LC,K1,2) 2909 end do 2910 ELSE 2911 DO K1=1,3 2912 DRCDC(K1,3)=DRDC(LC,K1,2) 2913 DRCDC(K1,4)=DRDC(LC,K1,1) 2914 end do 2915 ENDIF 2916 2917 do k1=1,3 2918 dhddc(1,k1)=dh1rdc(1,k1) 2919 dhddc(2,k1)=dh1rdc(2,k1) 2920 dhddc(3,k1)=dh1rdc(3,k1) 2921 dhedc(1,k1+1)=dh2rdc(1,k1) 2922 dhedc(2,k1+1)=dh2rdc(2,k1) 2923 dhedc(3,k1+1)=dh2rdc(3,k1) 2924 end do 2925 2926********************************************************************** 2927* write (64,*)j(1),j(2),j(3),j(4) 2928* do k1=1,3 2929* write (64,'(10f12.4)')(dh1rdc(k1,k2),k2=1,3), 2930* $(dhdc(ld,k1,k2),k2=1,3),(dhddc(k1,k2),k2=1,4) 2931* write (64,'(10f12.4)')(dh2rdc(k1,k2),k2=1,3), 2932* $(dhdc(le,k1,k2),k2=1,3),(dhedc(k1,k2),k2=1,4) 2933* end do 2934* write (64,*) 2935********************************************************************** 2936 HTRA=RLA+COSHD*(RLC*COSHE-RLB) 2937 HTRB=RLB-RLA*COSHD-RLC*COSHE 2938 HTRC=RLC+COSHE*(RLA*COSHD-RLB) 2939 HTHD=RLA*SINHD*(RLB-RLC*COSHE) 2940 HTHE=RLC*SINHE*(RLB-RLA*COSHD) 2941 HNRA=RLC*SINHD*SINHE 2942 HNRC=RLA*SINHD*SINHE 2943 HNHD=RLA*RLC*COSHD*SINHE 2944 HNHE=RLA*RLC*SINHD*COSHE 2945 2946 if (Lvirial.eq.1) then 2947 do k1=1,3 2948 do k2=1,3 2949 virial_tmp(k1,k2) = 0.0 2950 end do 2951 end do 2952 endif 2953 2954 DO K1=1,3 2955 DRDA(K1)=A(K1)/R4 2956 DO K2=1,4 2957 DRVDC(K1,K2)=DRDA(K1)*DADC(K2) 2958 DTDC(K1,K2)=2.0*(DRADC(K1,K2)*HTRA+DRBDC(K1,K2)*HTRB+DRCDC(K1,K2 2959 $)*HTRC-DRVDC(K1,K2)*R4+DHDDC(K1,K2)*HTHD+DHEDC(K1,K2)*HTHE) 2960 DNDC(K1,K2)=2.0*(DRADC(K1,K2)*HNRA+DRCDC(K1,K2)*HNRC+DHDDC(K1,K2 2961 $)*HNHD+DHEDC(K1,K2)*HNHE) 2962 DARGTDC(i1,K1,K2)=(DTDC(K1,K2)-ARG*DNDC(K1,K2))/POEM 2963 2964 ftmp = DARGTDC(i1,K1,K2)*detdar+ 2965 $dargtdc(i1,k1,k2)*decodar+(detdhd+decodhd)*dhddc(k1,k2)+ 2966 $(detdhe+decodhe)*dhedc(k1,k2) 2967 D(K1,J(K2))=D(K1,J(K2))+ftmp 2968 2969 if (Lvirial.eq.1) then 2970 do k1p=1,3 2971 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(j(k2),k1p) 2972 end do 2973 endif 2974 2975 end do 2976 end do 2977 if (Lvirial.eq.1) then 2978 virialsym(1) = virial_tmp(1,1) 2979 virialsym(2) = virial_tmp(2,2) 2980 virialsym(3) = virial_tmp(3,3) 2981 virialsym(4) = virial_tmp(1,2) 2982 virialsym(5) = virial_tmp(1,3) 2983 virialsym(6) = virial_tmp(2,3) 2984 do k1 = 1,6 2985 virial(k1) = virial(k1) + virialsym(k1) 2986 end do 2987 2988 if (Latomvirial.eq.1) then 2989 frac = 1.0d0/4 2990 do k1 = 1,6 2991 vtmp = virialsym(k1)*frac 2992 do k2=1,4 2993 ihu=j(k2) 2994 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 2995 end do 2996 end do 2997 endif 2998 2999 endif 3000 3001 if (Lvirial.eq.1) then 3002 do k1=1,3 3003 do k2=1,3 3004 virial_tmp(k1,k2) = 0.0 3005 end do 3006 end do 3007 endif 3008 3009 do i2=1,idbo1(la) 3010 ihu=idbo(la,i2) 3011 do k1=1,3 3012 ftmp = dbondc(la,k1,i2)*(detdboa+decodbola) 3013 d(k1,ihu)=d(k1,ihu)+ftmp 3014 3015 if (Lvirial.eq.1) then 3016 do k1p=1,3 3017 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 3018 end do 3019 endif 3020 3021 end do 3022 end do 3023 if (Lvirial.eq.1) then 3024 virialsym(1) = virial_tmp(1,1) 3025 virialsym(2) = virial_tmp(2,2) 3026 virialsym(3) = virial_tmp(3,3) 3027 virialsym(4) = virial_tmp(1,2) 3028 virialsym(5) = virial_tmp(1,3) 3029 virialsym(6) = virial_tmp(2,3) 3030 do k1 = 1,6 3031 virial(k1) = virial(k1) + virialsym(k1) 3032 end do 3033 3034 if (Latomvirial.eq.1) then 3035 frac = 1.0d0/idbo1(la) 3036 do k1 = 1,6 3037 vtmp = virialsym(k1)*frac 3038 do i2=1,idbo1(la) 3039 ihu=idbo(la,i2) 3040 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 3041 end do 3042 end do 3043 endif 3044 3045 endif 3046 3047 if (Lvirial.eq.1) then 3048 do k1=1,3 3049 do k2=1,3 3050 virial_tmp(k1,k2) = 0.0 3051 end do 3052 end do 3053 endif 3054 3055 do i2=1,idbo1(lb) 3056 ihu=idbo(lb,i2) 3057 do k1=1,3 3058 ftmp = dbondc(lb,k1,i2)*(detdbob+decodbolb) 3059 $ +dbopindc(lb,k1,i2)*detdbopib 3060 d(k1,ihu)=d(k1,ihu)+ftmp 3061 3062 if (Lvirial.eq.1) then 3063 do k1p=1,3 3064 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 3065 end do 3066 endif 3067 3068 end do 3069 end do 3070 if (Lvirial.eq.1) then 3071 virialsym(1) = virial_tmp(1,1) 3072 virialsym(2) = virial_tmp(2,2) 3073 virialsym(3) = virial_tmp(3,3) 3074 virialsym(4) = virial_tmp(1,2) 3075 virialsym(5) = virial_tmp(1,3) 3076 virialsym(6) = virial_tmp(2,3) 3077 do k1 = 1,6 3078 virial(k1) = virial(k1) + virialsym(k1) 3079 end do 3080 3081 if (Latomvirial.eq.1) then 3082 frac = 1.0d0/idbo1(lb) 3083 do k1 = 1,6 3084 vtmp = virialsym(k1)*frac 3085 do i2=1,idbo1(lb) 3086 ihu=idbo(lb,i2) 3087 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 3088 end do 3089 end do 3090 endif 3091 3092 endif 3093 3094 if (Lvirial.eq.1) then 3095 do k1=1,3 3096 do k2=1,3 3097 virial_tmp(k1,k2) = 0.0 3098 end do 3099 end do 3100 endif 3101 3102 do i2=1,idbo1(lc) 3103 ihu=idbo(lc,i2) 3104 do k1=1,3 3105 ftmp = dbondc(lc,k1,i2)*(detdboc+decodbolc) 3106 d(k1,ihu)=d(k1,ihu)+ftmp 3107 3108 if (Lvirial.eq.1) then 3109 do k1p=1,3 3110 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 3111 end do 3112 endif 3113 3114 end do 3115 end do 3116 if (Lvirial.eq.1) then 3117 virialsym(1) = virial_tmp(1,1) 3118 virialsym(2) = virial_tmp(2,2) 3119 virialsym(3) = virial_tmp(3,3) 3120 virialsym(4) = virial_tmp(1,2) 3121 virialsym(5) = virial_tmp(1,3) 3122 virialsym(6) = virial_tmp(2,3) 3123 do k1 = 1,6 3124 virial(k1) = virial(k1) + virialsym(k1) 3125 end do 3126 3127 if (Latomvirial.eq.1) then 3128 frac = 1.0d0/idbo1(lc) 3129 do k1 = 1,6 3130 vtmp = virialsym(k1)*frac 3131 do i2=1,idbo1(lc) 3132 ihu=idbo(lc,i2) 3133 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 3134 end do 3135 end do 3136 endif 3137 3138 endif 3139 3140 do i2=1,ia(j(2),2) 3141 iob=ia(j(2),2+i2) 3142 ncubo=nubon2(j(2),i2) 3143 if (bo(ncubo).gt.0.0) then 3144 3145 if (Lvirial.eq.1) then 3146 do k1=1,3 3147 do k2=1,3 3148 virial_tmp(k1,k2) = 0.0 3149 end do 3150 end do 3151 endif 3152 3153 do i3=1,idbo1(ncubo) 3154 ihu=idbo(ncubo,i3) 3155 do k1=1,3 3156 ftmp = detdsbo1*dbondc(ncubo,k1,i3) 3157 d(k1,ihu)=d(k1,ihu)+ftmp 3158 3159 if (Lvirial.eq.1) then 3160 do k1p=1,3 3161 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 3162 end do 3163 endif 3164 3165 end do 3166 end do 3167 if (Lvirial.eq.1) then 3168 virialsym(1) = virial_tmp(1,1) 3169 virialsym(2) = virial_tmp(2,2) 3170 virialsym(3) = virial_tmp(3,3) 3171 virialsym(4) = virial_tmp(1,2) 3172 virialsym(5) = virial_tmp(1,3) 3173 virialsym(6) = virial_tmp(2,3) 3174 do k1 = 1,6 3175 virial(k1) = virial(k1) + virialsym(k1) 3176 end do 3177 3178 if (Latomvirial.eq.1) then 3179 frac = 1.0d0/idbo1(ncubo) 3180 do k1 = 1,6 3181 vtmp = virialsym(k1)*frac 3182 do i3=1,idbo1(ncubo) 3183 ihu=idbo(ncubo,i3) 3184 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 3185 end do 3186 end do 3187 endif 3188 3189 endif 3190 3191 endif 3192 end do 3193 3194 do i2=1,ia(j(3),2) 3195 iob=ia(j(3),2+i2) 3196 ncubo=nubon2(j(3),i2) 3197 if (bo(ncubo).gt.0.0) then 3198 3199 if (Lvirial.eq.1) then 3200 do k1=1,3 3201 do k2=1,3 3202 virial_tmp(k1,k2) = 0.0 3203 end do 3204 end do 3205 endif 3206 3207 do i3=1,idbo1(ncubo) 3208 ihu=idbo(ncubo,i3) 3209 do k1=1,3 3210 ftmp = detdsbo1*dbondc(ncubo,k1,i3) 3211 d(k1,ihu)=d(k1,ihu)+ftmp 3212 3213 if (Lvirial.eq.1) then 3214 do k1p=1,3 3215 virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p) 3216 end do 3217 endif 3218 3219 end do 3220 end do 3221 3222 if (Lvirial.eq.1) then 3223 virialsym(1) = virial_tmp(1,1) 3224 virialsym(2) = virial_tmp(2,2) 3225 virialsym(3) = virial_tmp(3,3) 3226 virialsym(4) = virial_tmp(1,2) 3227 virialsym(5) = virial_tmp(1,3) 3228 virialsym(6) = virial_tmp(2,3) 3229 do k1 = 1,6 3230 virial(k1) = virial(k1) + virialsym(k1) 3231 end do 3232 3233 if (Latomvirial.eq.1) then 3234 frac = 1.0d0/idbo1(ncubo) 3235 do k1 = 1,6 3236 vtmp = virialsym(k1)*frac 3237 do i3=1,idbo1(ncubo) 3238 ihu=idbo(ncubo,i3) 3239 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 3240 end do 3241 end do 3242 endif 3243 3244 endif 3245 3246 endif 3247 end do 3248 3249 10 continue 3250 3251 return 3252 end 3253********************************************************************** 3254********************************************************************** 3255 3256 subroutine nonbon 3257 3258********************************************************************** 3259#include "cbka.blk" 3260#include "cbkc.blk" 3261#include "cbkch.blk" 3262#include "cbkconst.blk" 3263#include "cbkd.blk" 3264#include "cbkdcell.blk" 3265#include "cbkenergies.blk" 3266#include "cbkff.blk" 3267#include "cbkia.blk" 3268#include "cbknonbon.blk" 3269#include "cbkpairs.blk" 3270#include "cbknvlown.blk" 3271#include "cellcoord.blk" 3272#include "control.blk" 3273#include "small.blk" 3274 3275 dimension a(3),da(6) 3276 dimension virial_tmp(3,3),virialsym(6) 3277********************************************************************** 3278* * 3279* Calculate vdWaals and Coulomb energies and derivatives * 3280* * 3281********************************************************************** 3282c$$$ if (ndebug.eq.1) then 3283c$$$C open (65,file='fort.65',status='unknown',access='append') 3284c$$$ write (65,*) 'In nonbon' 3285c$$$ call timer(65) 3286c$$$ end if 3287 3288 ew=0.0 3289 ep=0.0 3290 3291 c1c=332.0638 3292 third=one/three 3293 fothird=4.0/3.0 3294 twothird=2.0/3.0 3295 h15=(vpar(29)-1.0)/vpar(29) 3296 3297 nptmp=0 3298 nstmp=0 3299 do 10 ivl=1,nvpair-nvlself 3300c Use precomputed midpoint criterion to decide if interaction is owned. 3301 if (nvlown(ivl).eq.1) then 3302 3303 i1=nvl1(ivl) 3304 i2=nvl2(ivl) 3305 3306 call dista2(i1,i2,rr,a(1),a(2),a(3)) 3307 if (rr.gt.swb.or.rr.lt.0.001) goto 10 3308 3309 ity1=ia(i1,1) 3310 ity2=ia(i2,1) 3311 imol1=iag(i1,3+mbond) 3312 imol2=iag(i2,3+mbond) 3313 rr2=rr*rr 3314 3315 sw=1.0 3316 sw1=0.0 3317 call taper(rr,rr2) 3318********************************************************************** 3319* * 3320* Calculate vdWaals energy * 3321* * 3322********************************************************************** 3323 p1=p1co(ity1,ity2) 3324 p2=p2co(ity1,ity2) 3325 p3=p3co(ity1,ity2) 3326 hulpw=(rr**vpar(29)+gamwco(ity1,ity2)) 3327 rrw=hulpw**(1.0/vpar(29)) 3328 h1=exp(p3*(1.0-rrw/p1)) 3329 h2=exp(0.50*p3*(1.0-rrw/p1)) 3330 3331 ewh=p2*(h1-2.0*h2) 3332 rrhuw=rr**(vpar(29)-1.0) 3333 dewdr=(p2*p3/p1)*(h2-h1)*rrhuw*(hulpw**(-h15)) 3334 3335********************************************************************** 3336* * 3337* Calculate Coulomb energy * 3338* * 3339********************************************************************** 3340 q1q2=ch(i1)*ch(i2) 3341 hulp1=(rr2*rr+gamcco(ity1,ity2)) 3342 eph=c1c*q1q2/(hulp1**third) 3343 depdr=-c1c*q1q2*rr2/(hulp1**fothird) 3344********************************************************************** 3345* * 3346* Taper correction * 3347* * 3348********************************************************************** 3349 ephtap=eph*sw 3350 depdrtap=depdr*sw+eph*sw1 3351 ewhtap=ewh*sw 3352 dewdrtap=dewdr*sw+ewh*sw1 3353 3354* write (64,*)i1,i2,p1,p2,p3,gamwco(ity1,ity2),vpar(29),ewh,ew 3355 ew=ew+ewhtap 3356 ep=ep+ephtap 3357 estrain(i1)=estrain(i1)+0.50*(ewhtap+ephtap) !1st atom energy 3358 estrain(i2)=estrain(i2)+0.50*(ewhtap+ephtap) !2nd atom energy 3359 3360********************************************************************** 3361* * 3362* Calculate derivatives vdWaals energy to cartesian * 3363* coordinates * 3364* * 3365********************************************************************** 3366 3367 if (Lvirial.eq.1) then 3368 do k1=1,3 3369 do k2=1,3 3370 virial_tmp(k1,k2) = 0.0 3371 end do 3372 end do 3373 endif 3374 3375 do k4=1,3 3376 ftmp = (dewdrtap+depdrtap)*(a(k4)/rr) 3377 d(k4,i1)=d(k4,i1)+ftmp 3378 d(k4,i2)=d(k4,i2)-ftmp 3379 if (Lvirial.eq.1) then 3380 do k1p=1,3 3381 virial_tmp(k4,k1p)=virial_tmp(k4,k1p)+ 3382 $ ftmp*c(i1,k1p)-ftmp*c(i2,k1p) 3383 end do 3384 endif 3385 end do 3386 3387 if (Lvirial.eq.1) then 3388 virialsym(1) = virial_tmp(1,1) 3389 virialsym(2) = virial_tmp(2,2) 3390 virialsym(3) = virial_tmp(3,3) 3391 virialsym(4) = virial_tmp(1,2) 3392 virialsym(5) = virial_tmp(1,3) 3393 virialsym(6) = virial_tmp(2,3) 3394 do k1 = 1,6 3395 virial(k1) = virial(k1) + virialsym(k1) 3396 end do 3397 3398 if (Latomvirial.eq.1) then 3399 frac = 1.0d0/2 3400 do k1 = 1,6 3401 vtmp = virialsym(k1)*frac 3402 ihu=i1 3403 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 3404 ihu=i2 3405 atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp 3406 end do 3407 endif 3408 3409 endif 3410 3411 endif 3412 3413 10 continue 3414 3415 return 3416 end 3417 3418********************************************************************** 3419********************************************************************** 3420 3421 subroutine efield 3422 3423********************************************************************** 3424#include "cbka.blk" 3425#include "cbkc.blk" 3426#include "cbkch.blk" 3427#include "cbkcha.blk" 3428#include "cbkd.blk" 3429#include "cbkefield.blk" 3430#include "cbkenergies.blk" 3431#include "cbktregime.blk" 3432#include "control.blk" 3433#include "small.blk" 3434c$$$ if (ndebug.eq.1) then 3435c$$$C open (65,file='fort.65',status='unknown',access='append') 3436c$$$ write (65,*) 'In efield' 3437c$$$ call timer(65) 3438c$$$ close (65) 3439c$$$ end if 3440********************************************************************** 3441* * 3442* Electric field * 3443* * 3444********************************************************************** 3445 efi=0.0 3446 efix=0.0 3447 efiy=0.0 3448 efiz=0.0 3449 c1c=332.0638 !Coulomb energy conversion 3450 c1=23.02 !conversion from kcal to eV 3451 3452 if (ifieldx.eq.1) then 3453 do i1=1,na 3454 efih=vfieldx*c1*c1c*ch(i1)*c(i1,1) 3455 efix=efix+efih 3456 estrain(i1)=estrain(i1)+efih !atom energy 3457 3458 defidc=c1*c1c*vfieldx*ch(i1) 3459 d(1,i1)=d(1,i1)+defidc 3460 end do 3461 end if 3462 3463 if (ifieldy.eq.1) then 3464 do i1=1,na 3465 efih=vfieldy*c1*c1c*ch(i1)*c(i1,2) 3466 efiy=efiy+efih 3467 estrain(i1)=estrain(i1)+efih !atom energy 3468 3469 defidc=c1*c1c*vfieldy*ch(i1) 3470 d(2,i1)=d(2,i1)+defidc 3471 end do 3472 end if 3473 3474 if (ifieldz.eq.1) then 3475 do i1=1,na 3476 efih=vfieldz*c1*c1c*ch(i1)*c(i1,3) 3477 efiz=efiz+efih 3478 estrain(i1)=estrain(i1)+efih !atom energy 3479 3480 defidc=c1*c1c*vfieldz*ch(i1) 3481 d(3,i1)=d(3,i1)+defidc 3482 end do 3483 end if 3484 3485 efi=efix+efiy+efiz 3486 return 3487 end 3488********************************************************************** 3489********************************************************************** 3490 3491 subroutine restraint 3492 3493********************************************************************** 3494#include "cbka.blk" 3495#include "cbkatomcoord.blk" 3496#include "cbkc.blk" 3497#include "cbkconst.blk" 3498#include "cbkd.blk" 3499#include "cbkenergies.blk" 3500#include "cbkrestr.blk" 3501#include "cbktorang.blk" 3502#include "cbktorsion.blk" 3503#include "cbktregime.blk" 3504#include "control.blk" 3505#include "small.blk" 3506#include "cbkinit.blk" 3507 dimension drda(3),j(4),dhrdc(3,3),dargdc(3,3) 3508********************************************************************** 3509* * 3510* Calculate distance restraint energy * 3511* * 3512********************************************************************** 3513c$$$ if (ndebug.eq.1) then 3514c$$$C open (65,file='fort.65',status='unknown',access='append') 3515c$$$ write (65,*) 'In restraint' 3516c$$$ call timer(65) 3517c$$$ close (65) 3518c$$$ end if 3519 do i1=1,nrestra 3520 ih1=irstra(i1,1) 3521 ih2=irstra(i1,2) 3522 if (itend(i1).eq.0.or.(mdstep.gt.itstart(i1).and.mdstep.lt. 3523 $itend(i1))) then 3524 call dista2(ih1,ih2,rr,dx,dy,dz) 3525 diffr=rr-rrstra(i1) 3526* diffr=rrstra(i1) 3527 exphu=exp(-vkrst2(i1)*(diffr*diffr)) 3528 erh=vkrstr(i1)*(1.0-exphu) 3529 deresdr=2.0*vkrst2(i1)*diffr*vkrstr(i1)*exphu 3530* deresdr=-2.0*vkrst2(i1)*diffr*vkrstr(i1)*exphu 3531 eres=eres+erh 3532 drda(1)=dx/rr 3533 drda(2)=dy/rr 3534 drda(3)=dz/rr 3535 do k1=1,3 3536 d(k1,ih1)=d(k1,ih1)+deresdr*drda(k1) 3537 d(k1,ih2)=d(k1,ih2)-deresdr*drda(k1) 3538 end do 3539 end if 3540 end do 3541 3542********************************************************************** 3543* * 3544* Calculate angle restraint energy * 3545* * 3546********************************************************************** 3547 do i1=1,nrestrav 3548 j(1)=irstrav(i1,1) 3549 j(2)=irstrav(i1,2) 3550 j(3)=irstrav(i1,3) 3551 ittr=0 3552* do i2=1,nval 3553* if (j(1).eq.iv(i2,2).and.j(2).eq.iv(i2,3).and.j(3).eq.iv(i2,4)) 3554* $ittr=i2 3555* end do 3556* if (ittr.eq.0) stop 'Wrong valence angle restraint' 3557 call calvalres(j(1),j(2),j(3),arg,hr,dhrdc,dargdc) 3558 vaval=hr*rdndgr 3559 diffv=-(vaval-vrstra(i1))*dgrrdn 3560 exphu=exp(-vkr2v(i1)*(diffv*diffv)) 3561 erh=vkrv(i1)*(1.0-exphu) 3562 deresdv=-2.0*vkr2v(i1)*diffv*vkrv(i1)*exphu 3563 eres=eres+erh 3564 do k1=1,3 3565 do k2=1,3 3566 d(k1,j(k2))=d(k1,j(k2))+deresdv*dhrdc(k1,k2) 3567 end do 3568 end do 3569 3570 end do 3571 3572********************************************************************** 3573* * 3574* Calculate torsion restraint energy * 3575* * 3576********************************************************************** 3577 do i1=1,nrestrat 3578 j(1)=irstrat(i1,1) 3579 j(2)=irstrat(i1,2) 3580 j(3)=irstrat(i1,3) 3581 j(4)=irstrat(i1,4) 3582 ittr=0 3583 do i2=1,ntor 3584 if (j(1).eq.it(i2,2).and.j(2).eq.it(i2,3).and.j(3).eq.it(i2,4) 3585 $.and.j(4).eq.it(i2,5)) ittr=i2 3586 if (j(4).eq.it(i2,2).and.j(3).eq.it(i2,3).and.j(2).eq.it(i2,4) 3587 $.and.j(1).eq.it(i2,5)) ittr=i2 3588 end do 3589 if (ittr.eq.0) then 3590 write (*,*)'Wrong torsion restraint' 3591 write (*,*)i1,j(1),j(2),j(3),j(4) 3592 stop 'Wrong torsion restraint' 3593 end if 3594 vtor=thg(ittr) 3595 difft=-(vtor-trstra(i1))*dgrrdn 3596 exphu=exp(-vkr2t(i1)*(difft*difft)) 3597 erh=vkrt(i1)*(1.0-exphu) 3598 deresdt=2.0*vkr2t(i1)*difft*vkrt(i1)*exphu 3599 if (vtor.lt.zero) deresdt=-deresdt 3600 eres=eres+erh 3601 do k1=1,3 3602 do k2=1,4 3603 d(k1,j(k2))=d(k1,j(k2))+deresdt*dargtdc(ittr,k1,k2) 3604 end do 3605 end do 3606 3607 end do 3608********************************************************************** 3609* * 3610* Calculate mass centre restraint energy * 3611* * 3612********************************************************************** 3613 do i1=1,nrestram 3614 j1=irstram(i1,2) 3615 j2=irstram(i1,3) 3616 j3=irstram(i1,4) 3617 j4=irstram(i1,5) 3618 kdir=irstram(i1,1) 3619 cmx1=0.0 3620 cmy1=0.0 3621 cmz1=0.0 3622 cmx2=0.0 3623 cmy2=0.0 3624 cmz2=0.0 3625 summas1=0.0 3626 summas2=0.0 3627 do i2=j1,j2 3628 cmx1=cmx1+c(i2,1)*xmasat(i2) 3629 cmy1=cmy1+c(i2,2)*xmasat(i2) 3630 cmz1=cmz1+c(i2,3)*xmasat(i2) 3631 summas1=summas1+xmasat(i2) 3632 end do 3633 cmx1=cmx1/summas1 3634 cmy1=cmy1/summas1 3635 cmz1=cmz1/summas1 3636 if (mdstep.lt.2) then 3637 rmstrax(i1)=cmx1 3638 rmstray(i1)=cmy1 3639 rmstraz(i1)=cmz1 3640 end if 3641 if (kdir.le.3) then 3642 do i2=j3,j4 3643 cmx2=cmx2+c(i2,1)*xmasat(i2) 3644 cmy2=cmy2+c(i2,2)*xmasat(i2) 3645 cmz2=cmz2+c(i2,3)*xmasat(i2) 3646 summas2=summas2+xmasat(i2) 3647 end do 3648 cmx2=cmx2/summas2 3649 cmy2=cmy2/summas2 3650 cmz2=cmz2/summas2 3651 end if 3652 if (kdir.eq.1) dist=cmx1-cmx2 3653 if (kdir.eq.2) dist=cmy1-cmy2 3654 if (kdir.eq.3) dist=cmz1-cmz2 3655 if (kdir.eq.4) then 3656 distx=cmx1-rmstrax(i1) 3657 disty=cmy1-rmstray(i1) 3658 distz=cmz1-rmstraz(i1) 3659 dist=sqrt(distx*distx+disty*disty+distz*distz) 3660 end if 3661 dismacen(i1)=dist 3662 dist=dist-rmstra1(i1) 3663 erh=rmstra2(i1)*dist*dist 3664 deresdr=2.0*dist*rmstra2(i1) 3665* exphu=exp(-rmstra3(i1)*(dist*dist)) 3666* erh=rmstra2(i1)*(1.0-exphu) 3667* deresdr=2.0*rmstra3(i1)*dist*rmstra2(i1)*exphu 3668 eres=eres+erh 3669 if (kdir.le.3) then 3670 do i2=j1,j2 3671 d(kdir,i2)=d(kdir,i2)+deresdr*xmasat(i2)/summas1 3672 end do 3673 do i2=j3,j4 3674 d(kdir,i2)=d(kdir,i2)-deresdr*xmasat(i2)/summas2 3675 end do 3676 end if 3677 if (kdir.eq.4.and.mdstep.gt.5) then 3678 do i2=j1,j2 3679 d(1,i2)=d(1,i2)+deresdr*(distx/dist)*(xmasat(i2)/summas1) 3680 d(2,i2)=d(2,i2)+deresdr*(disty/dist)*(xmasat(i2)/summas1) 3681 d(3,i2)=d(3,i2)+deresdr*(distz/dist)*(xmasat(i2)/summas1) 3682 end do 3683 end if 3684 end do 3685********************************************************************** 3686* * 3687* Calculate morphing energy * 3688* * 3689********************************************************************** 3690 if (imorph.eq.1) then 3691 distot=zero 3692 do i1=1,na 3693 dmx=c(i1,1)-cmo(i1,1) 3694 dmy=c(i1,2)-cmo(i1,2) 3695 dmz=c(i1,3)-cmo(i1,3) 3696 dism=sqrt(dmx*dmx+dmy*dmy+dmz*dmz) 3697 distot=distot+dism 3698* exphu=exp(-vmo2(i1)*(dism*dism)) 3699* erh=vmo1(i1)*(1.0-exphu) 3700 erh=vmo1(i1)*dism 3701 eres=eres+erh 3702* deresddis=2.0*vmo2(i1)*dism*vmo1(i1)*exphu 3703 deresddis=vmo1(i1) 3704 drda1=dmx/dism 3705 drda2=dmy/dism 3706 drda3=dmz/dism 3707 d(1,i1)=d(1,i1)+deresddis*drda1 3708 d(2,i1)=d(2,i1)+deresddis*drda2 3709 d(3,i1)=d(3,i1)+deresddis*drda3 3710 end do 3711 3712 end if 3713 3714 3715 return 3716 end 3717********************************************************************** 3718******************************************************************** 3719 3720 subroutine calvalres (ja1,ja2,ja3,arg,hr,dhrdc,dargdc) 3721 3722********************************************************************** 3723#include "cbka.blk" 3724#include "cbkc.blk" 3725 dimension a(3),b(3),j(3),dradc(3,3),drbdc(3,3),dtdc(3,3), 3726 $dargdc(3,3),dndc(3,3),dadc(3),dbdc(3),dhrdc(3,3) 3727********************************************************************** 3728* * 3729* Calculate valency angles and their derivatives to cartesian * 3730* coordinates for restraint calculations * 3731* * 3732********************************************************************** 3733c$$$* if (ndebug.eq.1) then 3734c$$$C* open (65,file='fort.65',status='unknown',access='append') 3735c$$$* write (65,*) 'In calvalres' 3736c$$$* call timer(65) 3737c$$$* close (65) 3738c$$$* end if 3739 3740 dadc(1)=-1.0 3741 dadc(2)=1.0 3742 dadc(3)=0.0 3743 dbdc(1)=0.0 3744 dbdc(2)=1.0 3745 dbdc(3)=-1.0 3746 do k1=1,3 3747 do k2=1,3 3748 dradc(k1,k2)=0.0 3749 drbdc(k1,k2)=0.0 3750 end do 3751 end do 3752********************************************************************** 3753* * 3754* Determine valency angle * 3755* * 3756********************************************************************** 3757 call dista2(ja1,ja2,rla,dx1,dy1,dz1) 3758 call dista2(ja2,ja3,rlb,dx2,dy2,dz2) 3759 3760 a(1)=-dx1 3761 a(2)=-dy1 3762 a(3)=-dz1 3763 b(1)=dx2 3764 b(2)=dy2 3765 b(3)=dz2 3766 poem=rla*rlb 3767 tel=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) 3768 arg=tel/poem 3769 arg2=arg*arg 3770 s1ma22=1.0-arg2 3771 if (s1ma22.lt.1.0d-10) s1ma22=1.0d-10 3772 s1ma2=sqrt(s1ma22) 3773 if (arg.gt.1.0) arg=1.0 3774 if (arg.lt.-1.0) arg=-1.0 3775 hr=acos(arg) 3776********************************************************************** 3777* * 3778* Calculate derivative valency angle to cartesian coordinates * 3779* * 3780********************************************************************** 3781 do k1=1,3 3782 dradc(k1,1)=-a(k1)/rla 3783 dradc(k1,2)=a(k1)/rla 3784 end do 3785 3786 do k1=1,3 3787 drbdc(k1,2)=b(k1)/rlb 3788 drbdc(k1,3)=-b(k1)/rlb 3789 end do 3790 3791 do k1=1,3 3792 do k2=1,3 3793 dndc(k1,k2)=rla*drbdc(k1,k2)+rlb*dradc(k1,k2) 3794 dtdc(k1,k2)=a(k1)*dbdc(k2)+b(k1)*dadc(k2) 3795 dargdc(k1,k2)=(dtdc(k1,k2)-arg*dndc(k1,k2))/poem 3796 dhrdc(k1,k2)=-dargdc(k1,k2)/s1ma2 3797 end do 3798 end do 3799 3800 10 continue 3801 3802 return 3803 end 3804********************************************************************** 3805******************************************************************** 3806 3807 subroutine calvalhb (ja1,ja2,ja3,ix,iy,iz,arg,hr,dhrdc,dargdc) 3808 3809********************************************************************** 3810#include "cbka.blk" 3811#include "cbkc.blk" 3812 dimension a(3),b(3),j(3),dradc(3,3),drbdc(3,3),dtdc(3,3), 3813 $dargdc(3,3),dndc(3,3),dadc(3),dbdc(3),dhrdc(3,3) 3814********************************************************************** 3815* * 3816* Calculate valency angles and their derivatives to cartesian * 3817* coordinates for hydrogen bond calculations * 3818* * 3819********************************************************************** 3820c$$$* if (ndebug.eq.1) then 3821c$$$* open (65,file='fort.65',status='unknown',access='append') 3822c$$$* write (65,*) 'In calvalhb' 3823c$$$* call timer(65) 3824c$$$* close (65) 3825c$$$* end if 3826 3827 dadc(1)=-1.0 3828 dadc(2)=1.0 3829 dadc(3)=0.0 3830 dbdc(1)=0.0 3831 dbdc(2)=1.0 3832 dbdc(3)=-1.0 3833 do k1=1,3 3834 do k2=1,3 3835 dradc(k1,k2)=0.0 3836 drbdc(k1,k2)=0.0 3837 end do 3838 end do 3839********************************************************************** 3840* * 3841* Determine valency angle * 3842* * 3843********************************************************************** 3844 call dista2(ja1,ja2,rla,dx1,dy1,dz1) 3845 call dista2(ja2,ja3,rlb,dx2,dy2,dz2) 3846 3847 a(1)=-dx1 3848 a(2)=-dy1 3849 a(3)=-dz1 3850 b(1)=dx2 3851 b(2)=dy2 3852 b(3)=dz2 3853 poem=rla*rlb 3854 tel=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) 3855 arg=tel/poem 3856 arg2=arg*arg 3857 s1ma22=1.0-arg2 3858 if (s1ma22.lt.1.0d-10) s1ma22=1.0d-10 3859 s1ma2=sqrt(s1ma22) 3860 if (arg.gt.1.0) arg=1.0 3861 if (arg.lt.-1.0) arg=-1.0 3862 hr=acos(arg) 3863********************************************************************** 3864* * 3865* Calculate derivative valency angle to cartesian coordinates * 3866* * 3867********************************************************************** 3868 do k1=1,3 3869 dradc(k1,1)=-a(k1)/rla 3870 dradc(k1,2)=a(k1)/rla 3871 end do 3872 3873 do k1=1,3 3874 drbdc(k1,2)=b(k1)/rlb 3875 drbdc(k1,3)=-b(k1)/rlb 3876 end do 3877 3878 do k1=1,3 3879 do k2=1,3 3880 dndc(k1,k2)=rla*drbdc(k1,k2)+rlb*dradc(k1,k2) 3881 dtdc(k1,k2)=a(k1)*dbdc(k2)+b(k1)*dadc(k2) 3882 dargdc(k1,k2)=(dtdc(k1,k2)-arg*dndc(k1,k2))/poem 3883 dhrdc(k1,k2)=-dargdc(k1,k2)/s1ma2 3884 end do 3885 end do 3886 3887 10 continue 3888 3889 return 3890 end 3891********************************************************************** 3892********************************************************************** 3893 3894 subroutine caltor(ja1,ja2,ja3,ja4,ht) 3895 3896********************************************************************** 3897#include "cbka.blk" 3898#include "cbkenergies.blk" 3899#include "cbktregime.blk" 3900#include "control.blk" 3901#include "cbkinit.blk" 3902 DIMENSION A(3),DRDA(3),DADC(4),DRADC(3,4),DRBDC(3,4), 3903 $DRCDC(3,4),DHDDC(3,4),DHEDC(3,4),DRVDC(3,4),DTDC(3,4), 3904 $DNDC(3,4) 3905 dimension j(4),dvdc1(3,3),dargdc1(3,3),dvdc2(3,3),dargdc2(3,3) 3906********************************************************************** 3907* * 3908* Calculate torsion angle (for internal coordinates output) * 3909* * 3910********************************************************************** 3911c$$$ if (ndebug.eq.1) then 3912c$$$C open (65,file='fort.65',status='unknown',access='append') 3913c$$$ write (65,*) 'In caltor' 3914c$$$ call timer(65) 3915c$$$ close (65) 3916c$$$ end if 3917 do k1=1,3 3918 do k2=1,4 3919 dhddc(k1,k2)=0.0 3920 dhedc(k1,k2)=0.0 3921 dradc(k1,k2)=0.0 3922 drbdc(k1,k2)=0.0 3923 drcdc(k1,k2)=0.0 3924 end do 3925 end do 3926 et=0.0 3927 eco=0.0 3928 dadc(1)=1.0 3929 dadc(2)=0.0 3930 dadc(3)=0.0 3931 dadc(4)=-1.0 3932 call dista2(ja1,ja2,rla,dx1,dy1,dz1) 3933 call dista2(ja2,ja3,rlb,dx2,dy2,dz2) 3934 call dista2(ja3,ja4,rlc,dx2,dy2,dz2) 3935 call dista2(ja1,ja4,r4,dx2,dy2,dz2) 3936 call calvalres(ja1,ja2,ja3,arg1,h1,dvdc1,dargdc1) 3937 call calvalres(ja2,ja3,ja4,arg2,h2,dvdc2,dargdc2) 3938********************************************************************** 3939* * 3940* Determine torsion angle * 3941* * 3942********************************************************************** 3943 d142=r4*r4 3944 coshd=cos(h1) 3945 coshe=cos(h2) 3946 sinhd=sin(h1) 3947 sinhe=sin(h2) 3948 poem=2.0*rla*rlc*sinhd*sinhe 3949 poem2=poem*poem 3950 tel=rla*rla+rlb*rlb+rlc*rlc-d142-2.0*(rla*rlb*coshd-rla*rlc* 3951 $coshd*coshe+rlb*rlc*coshe) 3952 arg=tel/poem 3953 if (arg.gt.1.0) arg=1.0 3954 if (arg.lt.-1.0) arg=-1.0 3955 arg2=arg*arg 3956 ht=acos(arg)*rdndgr 3957 k1=ja1 3958 k2=ja2 3959 k3=ja3 3960 k4=ja4 3961 call dista2(k3,k2,dis,x3,y3,z3) 3962 y32z32=y3*y3+z3*z3 3963 wort1=sqrt(y32z32)+1e-6 3964 wort2=sqrt(y32z32+x3*x3)+1e-6 3965 sinalf=y3/wort1 3966 cosalf=z3/wort1 3967 sinbet=x3/wort2 3968 cosbet=wort1/wort2 3969 call dista2(k1,k2,dis,x1,y1,z1) 3970 x1=x1*cosbet-y1*sinalf*sinbet-z1*cosalf*sinbet 3971 y1=y1*cosalf-z1*sinalf 3972 wort3=sqrt(x1*x1+y1*y1)+1e-6 3973 singam=y1/wort3 3974 cosgam=x1/wort3 3975 call dista2(k4,k2,dis,x4,y4,z4) 3976 x4=x4*cosbet-y4*sinalf*sinbet-z4*cosalf*sinbet 3977 y4=y4*cosalf-z4*sinalf 3978 y4=x4*singam-y4*cosgam 3979 if (y4.gt.0.0) ht=-ht 3980 if (ht.lt.-179.999999d0) ht=-179.999999d0 3981 if (ht.gt.179.999999d0) ht=179.999999d0 3982 3983 return 3984 end 3985********************************************************************** 3986