1! --- 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt . 6! See Docs/Contributors.txt for a list of contributors. 7! --- 8 module molecularmechanics 9! 10! Add additional interactions through pair-potentials 11! Implementation: Julian Gale (Curtin, AU) 12! Modified by Alberto Garcia (stress sign fix, plots of V(r)) 13! Implemented Grimme's cutoff function (A. Garcia, April 2008) 14! Updated for arbitrary number of potentials (N. Papior, 2016) 15! 16 use precision, only: dp 17 implicit none 18 19! The original implementation by Julian had the wrong sign 20! for the stress tensor. 21! Define sign_change as +1 to recover that behavior 22 integer, parameter :: sign_change = -1 23 24 private 25 26 integer, save :: nMMpot = 0 27 integer, pointer, save :: MMpotptr(:,:) => null() ! (2,:) 28 integer, pointer, save :: MMpottype(:) => null() ! (:) 29 real(dp), save :: MMcutoff 30 real(dp), save :: s6_grimme 31 real(dp), save :: d_grimme 32 real(dp), pointer, save :: MMpotpar(:,:) => null() ! (6,:) 33 34 public :: inittwobody, twobody, reset_twobody 35 36 CONTAINS 37 38 subroutine inittwobody() 39 use fdf 40 use sys, only : die 41 use parallel,only : IONode 42 43 character(len=80) :: potential 44 character(len=80) :: scale 45 integer :: sizeMMpot 46 integer :: ni 47 integer :: nn 48 integer :: nr 49 real(dp) :: Lscale 50 real(dp) :: Escale 51 52 type(block_fdf) :: bfdf 53 type(parsed_line), pointer :: pline 54 55 ! initialize all the arrays, and allocate them dynamically. 56 call prep_arrays() 57 58 59! Get potential cutoff 60 MMcutoff = fdf_get('MM.Cutoff',30._dp,'Bohr') 61 62! Set MM units of energy for potential parameters 63 scale = fdf_get('MM.UnitsEnergy', 'eV') 64 ! this will allow arbitrary energy conversions 65 Escale = fdf_convfac(scale, 'Ry') 66 scale = fdf_get('MM.UnitsDistance', 'Ang') 67 ! this will allow arbitrary length conversions 68 Lscale = fdf_convfac(scale, 'Bohr') 69 70 ! If there are no molecular potentials we may quickly 71 ! exit 72 if ( sizeMMpot == 0 ) then 73 nMMpot = 0 74 return 75 end if 76 77 78! Read in data from block 79 nMMpot = 0 ! reset to count the actual number of potentials 80 if ( fdf_block('MM.Potentials', bfdf) ) then 81 if (IONode) then 82 write(6,"(a)") "Reading two-body potentials" 83 endif 84 85! Read and parse data line 86 do while ( fdf_bline(bfdf,pline) ) 87 88 nn = fdf_bnnames(pline) 89 ! If there are no names on this line 90 ! there is no specification of the type of 91 ! potential used. Hence, we immediately skip it 92 if ( nn == 0 ) cycle 93 94 ni = fdf_bnintegers(pline) 95 nr = fdf_bnreals(pline) 96 potential = fdf_bnames(pline,1) 97 98 if (leqi(potential,'C6')) then 99 call step_species(nMMpot, 1, ni, 'C6 - two-body') 100 101 if (nr .ge. 2) then 102! C6 : Parameter one is C6 coefficient 103 MMpotpar(1,nMMpot) = fdf_breals(pline,1)*Escale*(Lscale**6) 104! C6 : Parameter two is damping exponent 105 MMpotpar(2,nMMpot) = fdf_breals(pline,2)/Lscale 106 else if (nr .eq. 1) then 107! C6 : Parameter one is C6 coefficient 108 MMpotpar(1,nMMpot) = fdf_breals(pline,1)*Escale*(Lscale**6) 109 MMpotpar(2,nMMpot) = 0.0_dp 110 end if 111 112 else if (leqi(potential,'C8')) then 113 call step_species(nMMpot, 2, ni, 'C8 - two-body') 114 115 if (nr .ge. 2) then 116! C8 : Parameter one is C8 coefficient 117 MMpotpar(1,nMMpot) = fdf_breals(pline,1)*Escale*(Lscale**6) 118! C8 : Parameter two is damping exponent 119 MMpotpar(2,nMMpot) = fdf_breals(pline,2)/Lscale 120 elseif (nr .eq. 1) then 121! C8 : Parameter one is C8 coefficient 122 MMpotpar(1,nMMpot) = fdf_breals(pline,1)*Escale*(Lscale**6) 123 MMpotpar(2,nMMpot) = 0.0_dp 124 end if 125 126 else if (leqi(potential,'C10')) then 127 call step_species(nMMpot, 3, ni, 'C10 - two-body') 128 129 if (nr .ge. 2) then 130! C10 : Parameter one is C10 coefficient 131 MMpotpar(1,nMMpot) = fdf_breals(pline,1)*Escale*(Lscale**6) 132! C10 : Parameter two is damping exponent 133 MMpotpar(2,nMMpot) = fdf_breals(pline,2)/Lscale 134 else if (nr.eq.1) then 135! C10 : Parameter one is C10 coefficient 136 MMpotpar(1,nMMpot) = fdf_breals(pline,1)*Escale*(Lscale**6) 137 MMpotpar(2,nMMpot) = 0.0_dp 138 end if 139 140 else if (leqi(potential,'HARM')) then 141 call step_species(nMMpot, 4, ni, 'Harmonic - two-body') 142 143 if (nr .ge. 2) then 144! Harm : Parameter one is force constant 145 MMpotpar(1,nMMpot) = fdf_breals(pline,1)*Escale/(Lscale**2) 146! Harm : Parameter two is r0 147 MMpotpar(2,nMMpot) = fdf_breals(pline,2)*Lscale 148 else if (nr .eq. 1) then 149! Harm : Parameter one is force constant 150 MMpotpar(1,nMMpot) = fdf_breals(pline,1)*Escale/(Lscale**2) 151 MMpotpar(2,nMMpot) = 0.0_dp 152 end if 153 154 else if (leqi(potential,'Grimme')) then 155 call step_species(nMMpot, 5, ni, 'Grimme - two-body') 156 157 if (nr.eq.2) then 158! C6 : Parameter one is C6 coefficient 159 MMpotpar(1,nMMpot) = fdf_breals(pline,1)*Escale*(Lscale**6) 160 161! C6 : Parameter two is the sum of the van-der-Waals radii 162! Note 1: This must be already appropriately corrected (i.e., factor of 1.1 ...) 163! Note 2: This is a real length, as opposed to the damping parameters for the Tang-Toenes 164! potentials, so note the correct application of the scale factor. 165 MMpotpar(2,nMMpot) = fdf_breals(pline,2) * Lscale 166 167 else 168 call die('MM: Need both C6 and R0 values in Grimme line!') 169 end if 170 end if 171 172 end do 173 174 ! R. Peverati for DZ basis sets 175 s6_grimme = fdf_get("MM.Grimme.S6",1.66_dp) 176 ! 2006 Grimme paper 177 d_grimme = fdf_get("MM.Grimme.D",20._dp) 178 179 if (IONode) call plot_functions() 180 181 endif 182 183 if ( sizeMMpot /= nMMpot ) then 184 call die("MM: Too many lines in MM.Potentials block are not & 185 &read. Please delete non-applicable lines.") 186 end if 187 188 contains 189 190 !> @param MMpot current potential index 191 !> @param method the method value for the current potential 192 !> @param ni number of integers on the line (simple check) 193 !> @param name additional print out of method name 194 subroutine step_species(MMpot, method, ni, name) 195 integer, intent(inout) :: MMpot 196 integer, intent(in) :: method, ni 197 character(len=*), intent(in) :: name 198 199 ! Step potential 200 MMpot = MMpot + 1 201 202 ! Assign potential type 203 MMpottype(MMpot) = method 204 205 ! Check whether the correct species are found 206 if ( ni < 2 ) then 207 call die('MM: Species numbers missing in potential input!') 208 end if 209 210 ! Read the species indices 211 MMpotptr(1,MMpot) = fdf_bintegers(pline,1) 212 MMpotptr(2,MMpot) = fdf_bintegers(pline,2) 213 if ( IONode ) then 214 write(*,"(a,2(a,i0))") name, " potential between ", & 215 MMpotptr(1,MMpot), " and ", MMpotptr(2,MMpot) 216 end if 217 218 end subroutine step_species 219 220 ! dynamically create the arrays to allow 221 ! arbitrary number of potentials attached. 222 subroutine prep_arrays() 223 use alloc, only : re_alloc 224 type(block_fdf) :: bfdf 225 type(parsed_line), pointer :: pline 226 integer :: nn,ni,nr 227 character(len=80) :: potential 228 229 ! the global number of potentials 230 sizeMMpot = 0 231 if ( .not. fdf_block('MM.Potentials',bfdf) ) return 232 do while ( fdf_bline(bfdf,pline) ) 233 234 ! ensure that we can read the content 235 nn = fdf_bnnames(pline) 236 potential = fdf_bnames(pline,1) 237 ! we have a potential line... 238 ! note that if it an invalid line, this will 239 ! assert that we have (allocated potentials) > (actual potentials) 240 if ( nn > 0 ) sizeMMpot = sizeMMpot + 1 241 242 end do 243 244 if ( sizeMMpot == 0 ) return 245 246 ! allocate them 247 call re_alloc(MMpotptr , 1 , 2 , 1 , sizeMMpot, & 248 'MMpotptr','twobody') 249 call re_alloc(MMpottype , 1 , sizeMMpot, & 250 'MMpottype','twobody') 251 call re_alloc(MMpotpar , 1 , 6 , 1 , sizeMMpot, & 252 'MMpotpar','twobody') 253 254 ! initialize them to zero 255 ! ensures that we don't accidentally assign wrong 256 ! potential to a specie 257 MMpotptr(:,:) = 0 258 MMpottype(:) = 0 259 MMpotpar(:,:) = 0._dp 260 261 end subroutine prep_arrays 262 263 end subroutine inittwobody 264 265 subroutine reset_twobody() 266 267 use alloc, only: de_alloc 268 269 ! deallocate them 270 call de_alloc(MMpotptr, 'MMpotptr','twobody') 271 call de_alloc(MMpottype, 'MMpottype','twobody') 272 call de_alloc(MMpotpar, 'MMpotpar','twobody') 273 nMMpot = 0 274 275 end subroutine reset_twobody 276 277 278 subroutine twobody( na, xa, isa, cell, emm, ifa, fa, istr, stress ) 279 use parallel, only : IONode 280 use units, only : kbar 281 use alloc, only : re_alloc, de_alloc 282 283 integer, intent(in) :: na 284 integer, intent(in) :: isa(*) 285 integer, intent(in) :: ifa ! compute forces if > 0 286 integer, intent(in) :: istr ! compute stress if > 0 287 real(dp), intent(in) :: cell(3,3) 288 real(dp), intent(in) :: xa(3,*) 289 real(dp), intent(out) :: emm 290 real(dp), intent(inout) :: fa(3,*) 291 real(dp), intent(inout) :: stress(3,3) 292 293 integer :: i 294 integer :: idir 295 integer :: ii 296 integer :: imid 297 integer :: j 298 integer :: jdir 299 integer :: jj 300 integer :: jmid 301 integer :: kdir 302 integer :: kk 303 integer :: kmid 304 integer :: np 305 integer :: nvec 306 integer :: nvecj0 307 integer :: nveck0 308 logical :: lallfound1 309 logical :: lallfound2 310 logical :: lallfound3 311 logical :: lanyvalidpot 312 logical, pointer :: lvalidpot(:) 313 real(dp) :: MMcutoff2 314 real(dp) :: arg 315 real(dp) :: a 316 real(dp) :: b 317 real(dp) :: br6 318 real(dp) :: br8 319 real(dp) :: br10 320 real(dp) :: c 321 real(dp) :: alpha 322 real(dp) :: beta 323 real(dp) :: df6 324 real(dp) :: df8 325 real(dp) :: df10 326 real(dp) :: gamma 327 real(dp) :: earg 328 real(dp) :: ebr6 329 real(dp) :: ebr8 330 real(dp) :: ebr10 331 real(dp) :: etrm 332 real(dp) :: etrm1 333 real(dp) :: f6 334 real(dp) :: f8 335 real(dp) :: f10 336 real(dp) :: fg 337 real(dp) :: fgprime 338 real(dp) :: ftrm 339 real(dp) :: factor 340 real(dp) :: proj1 341 real(dp) :: proj2 342 real(dp) :: proj3 343 real(dp) :: r 344 real(dp) :: R0 345 real(dp) :: r2 346 real(dp) :: r2i 347 real(dp) :: r2j 348 real(dp) :: r2k 349 real(dp) :: rcx1 350 real(dp) :: rcy1 351 real(dp) :: rcz1 352 real(dp) :: rcx2 353 real(dp) :: rcy2 354 real(dp) :: rcz2 355 real(dp) :: rcx3 356 real(dp) :: rcy3 357 real(dp) :: rcz3 358 real(dp) :: recipa 359 real(dp) :: recipb 360 real(dp) :: recipc 361 real(dp) :: rnorm 362 real(dp) :: rx 363 real(dp) :: ry 364 real(dp) :: rz 365 real(dp) :: rxi 366 real(dp) :: ryi 367 real(dp) :: rzi 368 real(dp) :: rxj 369 real(dp) :: ryj 370 real(dp) :: rzj 371 real(dp) :: rvol 372 real(dp) :: vol 373 real(dp), external :: volcel 374 real(dp) :: x 375 real(dp) :: y 376 real(dp) :: z 377 378 real(dp) :: mm_stress(3,3) 379 integer :: jx 380 381! Start timer 382 call timer('MolMec', 1 ) 383 384! Initialise energy and mm_stress 385 emm = 0.0_dp 386 mm_stress(1:3,1:3) = 0.0_dp 387 388! Find number of cell images required 389 MMcutoff2 = MMcutoff**2 390 call uncell(cell,a,b,c,alpha,beta,gamma,1.0_dp) 391 recipa = 1.0_dp/a 392 recipb = 1.0_dp/b 393 recipc = 1.0_dp/c 394 395! Find volume if required 396 if (istr.ne.0) then 397 vol = volcel(cell) 398 rvol = 1.0_dp/vol 399 endif 400 401 ! quick return if no potentials are present 402 if ( nMMpot == 0 ) then 403 call timer('MolMec',2) 404 return 405 end if 406 407! Allocate workspace arrays 408 nullify(lvalidpot) 409 call re_alloc( lvalidpot, 1, nMMpot, 'lvalidpot', 'twobody' ) 410 411! Loop over first atom 412 do i = 1,na 413! Loop over second atom 414 do j = 1,i 415 if (i.eq.j) then 416 factor = 0.5_dp 417 else 418 factor = 1.0_dp 419 endif 420 421! Find valid potentials 422 lanyvalidpot = .false. 423 do np = 1,nMMpot 424 if ( MMpotptr(1,np) == isa(i) .and. & 425 MMpotptr(2,np) == isa(j)) then 426 lanyvalidpot = .true. 427 lvalidpot(np) = .true. 428 elseif ( MMpotptr(1,np) == isa(j) .and. & 429 MMpotptr(2,np) == isa(i)) then 430 lanyvalidpot = .true. 431 lvalidpot(np) = .true. 432 else 433 lvalidpot(np) = .false. 434 endif 435 enddo 436 if (lanyvalidpot) then 437! Find image of j nearest to i 438 x = xa(1,j) - xa(1,i) 439 y = xa(2,j) - xa(2,i) 440 z = xa(3,j) - xa(3,i) 441 442! Find projection of cell vector 3 on to i - j vector 443 rnorm = x*x + y*y + z*z 444 if (rnorm.gt.1.0d-12) rnorm = 1.0_dp/sqrt(rnorm) 445 proj3 = rnorm*recipc*(x*cell(1,3) + y*cell(2,3) + z*cell(3,3)) 446 kmid = nint(proj3) 447 x = x - kmid*cell(1,3) 448 y = y - kmid*cell(2,3) 449 z = z - kmid*cell(3,3) 450 451! Find projection of cell vector 2 on to i - j vector 452 rnorm = x*x + y*y + z*z 453 if (rnorm.gt.1.0d-12) rnorm = 1.0_dp/sqrt(rnorm) 454 proj2 = rnorm*recipb*(x*cell(1,2) + y*cell(2,2) + z*cell(3,2)) 455 jmid = nint(proj2) 456 x = x - jmid*cell(1,2) 457 y = y - jmid*cell(2,2) 458 z = z - jmid*cell(3,2) 459 460! Find projection of cell vector 1 on to i - j vector 461 rnorm = x*x + y*y + z*z 462 if (rnorm.gt.1.0d-12) rnorm = 1.0_dp/sqrt(rnorm) 463 proj1 = rnorm*recipa*(x*cell(1,1) + y*cell(2,1) + z*cell(3,1)) 464 imid = nint(proj1) 465 x = x - imid*cell(1,1) 466 y = y - imid*cell(2,1) 467 z = z - imid*cell(3,1) 468 469! Initialise counter for number of valid vectors 470 nvec = 0 471 472! Outer loop over first cell vector direction 473 do idir = 1,-1,-2 474! Reinitialise distance squared 475 r2i = 10000.0_dp*MMcutoff2 476 477! Loop over first cell vector 478 lallfound1 = .false. 479 if (idir == 1) then 480 ii = 0 481 else 482 ii = - 1 483 endif 484 485! Set initial coordinate vector 486 rxi = x + dble(ii)*cell(1,1) 487 ryi = y + dble(ii)*cell(2,1) 488 rzi = z + dble(ii)*cell(3,1) 489 490! Set increment vector 491 rcx1 = dble(idir)*cell(1,1) 492 rcy1 = dble(idir)*cell(2,1) 493 rcz1 = dble(idir)*cell(3,1) 494 495 do while (.not.lallfound1) 496! Save number of vectors before search over second direction 497 nvecj0 = nvec 498 499! Outer loop over second cell vector direction 500 do jdir = 1,-1,-2 501! Reinitialise saved distance squared 502 r2j = 10000.0_dp*MMcutoff2 503 504! Loop over second cell vector 505 lallfound2 = .false. 506 if (jdir == 1) then 507 jj = 0 508 else 509 jj = - 1 510 endif 511 512! Set initial coordinate vector 513 rxj = rxi + dble(jj)*cell(1,2) 514 ryj = ryi + dble(jj)*cell(2,2) 515 rzj = rzi + dble(jj)*cell(3,2) 516 517! Set increment vector 518 rcx2 = dble(jdir)*cell(1,2) 519 rcy2 = dble(jdir)*cell(2,2) 520 rcz2 = dble(jdir)*cell(3,2) 521 522 do while (.not.lallfound2) 523! Save number of vectors before search over third direction 524 nveck0 = nvec 525 526! Outer loop over third cell vector direction 527 do kdir = 1,-1,-2 528! Reinitialise saved distance squared 529 r2k = 10000.0_dp*MMcutoff2 530 531! Loop over third cell vector 532 lallfound3 = .false. 533 if (kdir == 1) then 534 kk = 0 535 else 536 kk = - 1 537 endif 538 539! Set initial coordinate vector 540 rx = rxj + dble(kk)*cell(1,3) 541 ry = ryj + dble(kk)*cell(2,3) 542 rz = rzj + dble(kk)*cell(3,3) 543 544! Set increment vector 545 rcx3 = dble(kdir)*cell(1,3) 546 rcy3 = dble(kdir)*cell(2,3) 547 rcz3 = dble(kdir)*cell(3,3) 548 549 do while (.not.lallfound3) 550! Calculate square of distance 551 r2 = rx*rx + ry*ry + rz*rz 552 553! Check distance squared against cutoff squared 554 if (r2.le.MMcutoff2) then 555 if (r2.gt.1.0d-10) then 556! Valid distance, so increment counter 557 nvec = nvec + 1 558! Evaluate potentials for this valid distance 559 do np = 1,nMMpot 560 if (lvalidpot(np)) then 561 select case ( MMpottype(np) ) 562 case ( 1 ) ! C6 563 if (MMpotpar(2,np) == 0.0_dp) then 564 etrm = MMpotpar(1,np)/(r2**3) 565 ftrm = 6.0_dp*factor*etrm/r2 566 etrm = - etrm 567 else 568 br6 = MMpotpar(2,np)*sqrt(r2) 569 ebr6 = exp(-br6) 570 f6 = 1.0_dp + br6*(1.0_dp + 0.5_dp*br6*(1.0_dp + (br6/3.0_dp)*( & 571 1.0_dp + 0.25_dp*br6*(1.0_dp + 0.2_dp*br6*(1.0_dp + (br6/6.0_dp)))))) 572 f6 = 1.0_dp - f6*ebr6 573 etrm1 = MMpotpar(1,np)/(r2**3) 574 df6 = ebr6*(br6*(br6**6))/720.0_dp 575 etrm = - etrm1*f6 576 ftrm = factor*(6.0_dp*etrm1*f6 - etrm1*df6)/r2 577 endif 578 case ( 2 ) ! C8 579 if (MMpotpar(2,np) == 0.0_dp) then 580 etrm = MMpotpar(1,np)/(r2**4) 581 ftrm = 8.0_dp*factor*etrm/r2 582 etrm = - etrm 583 else 584 br8 = MMpotpar(2,np)*sqrt(r2) 585 ebr8 = exp(-br8) 586 f8 = 1.0_dp + br8*(1.0_dp + 0.5_dp*br8*(1.0_dp + (br8/3.0_dp)*( & 587 1.0_dp + 0.25_dp*br8*(1.0_dp + 0.2_dp*br8*(1.0_dp + & 588 (br8/6.0_dp)*(1.0_dp+(br8/7.0_dp)*(1.0_dp + 0.125_dp*br8))))))) 589 f8 = 1.0_dp - f8*ebr8 590 etrm1 = MMpotpar(1,np)/(r2**4) 591 df8 = ebr8*(br8*(br8**8))/40320.0_dp 592 etrm = - etrm1*f8 593 ftrm = factor*(8.0_dp*etrm1*f8 - etrm1*df8)/r2 594 endif 595 case ( 3 ) ! C10 596 if (MMpotpar(2,np) == 0.0_dp) then 597 etrm = MMpotpar(1,np)/(r2**5) 598 ftrm = 10.0_dp*factor*etrm/r2 599 etrm = - etrm 600 else 601 br10 = MMpotpar(2,np)*sqrt(r2) 602 ebr10 = exp(-br10) 603 f10 = 1.0_dp + br10*(1.0_dp + 0.5_dp*br10*(1.0_dp + & 604 (br10/3.0_dp)*(1.0_dp + 0.25_dp*br10*(1.0_dp + 0.2_dp*br10*( & 605 1.0_dp + (br10/6.0_dp)*(1.0_dp + (br10/7.0_dp)*(1.0_dp + & 606 0.125_dp*br10*(1.0_dp + (br10/9.0_dp)*(1.0_dp + 0.1_dp*br10))))))))) 607 f10 = 1.0_dp - f10*ebr10 608 etrm1 = MMpotpar(1,np)/(r2**5) 609 df10 = ebr10*(br10*(br10**10))/3628800.0_dp 610 etrm = - etrm1*f10 611 ftrm = factor*(10.0_dp*etrm1*f10 - etrm1*df10)/r2 612 endif 613 case ( 4 ) ! Harm 614 r = sqrt(r2) 615 etrm = MMpotpar(1,np)*(r - MMpotpar(2,np)) 616 ftrm = factor*etrm/r 617 etrm = 0.5_dp*etrm*(r - MMpotpar(2,np)) 618 case ( 5 ) ! Grimme 619 r = sqrt(r2) 620 R0 = MMpotpar(2,np) 621 arg = - d_grimme*(r/R0 - 1.0_dp) 622 earg = exp(arg) 623 fg = 1.0_dp / ( 1.0_dp + earg ) 624 etrm1 = s6_grimme * MMpotpar(1,np)/(r2**3) 625 fgprime = (d_grimme/R0) * earg / ( 1 + earg )**2 626 etrm = - etrm1*fg 627 ftrm = factor*(6.0_dp*etrm1*fg - etrm1*r*fgprime)/r2 628 end select 629 630 emm = emm + factor*etrm 631 if (ifa.ne.0) then 632! Gradients 633 fa(1,i) = fa(1,i) + ftrm*rx 634 fa(2,i) = fa(2,i) + ftrm*ry 635 fa(3,i) = fa(3,i) + ftrm*rz 636 fa(1,j) = fa(1,j) - ftrm*rx 637 fa(2,j) = fa(2,j) - ftrm*ry 638 fa(3,j) = fa(3,j) - ftrm*rz 639 endif 640 if (istr.ne.0) then 641! Stress 642 ftrm = ftrm*rvol 643! Sign_change should be -1 for the correct 644! stress. See parameter definition above 645 mm_stress(1,1) = mm_stress(1,1) - & 646 sign_change*ftrm*rx*rx 647 mm_stress(2,1) = mm_stress(2,1) - & 648 sign_change*ftrm*ry*rx 649 mm_stress(3,1) = mm_stress(3,1) - & 650 sign_change*ftrm*rz*rx 651 mm_stress(1,2) = mm_stress(1,2) - & 652 sign_change*ftrm*rx*ry 653 mm_stress(2,2) = mm_stress(2,2) - & 654 sign_change*ftrm*ry*ry 655 mm_stress(3,2) = mm_stress(3,2) - & 656 sign_change*ftrm*rz*ry 657 mm_stress(1,3) = mm_stress(1,3) - & 658 sign_change*ftrm*rx*rz 659 mm_stress(2,3) = mm_stress(2,3) - & 660 sign_change*ftrm*ry*rz 661 mm_stress(3,3) = mm_stress(3,3) - & 662 sign_change*ftrm*rz*rz 663 endif ! istr.ne.0 664 endif ! lvalidpot(np) 665 enddo ! End loop nMMpot 666 endif ! Endif .gt. 1.0d-10 667 endif ! Endif .lt. MMcutoff2 668! Increment by third vector 669 kk = kk + kdir 670 rx = rx + rcx3 671 ry = ry + rcy3 672 rz = rz + rcz3 673 674! Check to see if this direction is complete 675 lallfound3 = (r2.gt.r2k.and.r2.gt.MMcutoff2) 676 r2k = r2 677 enddo ! End loop lallfound3 678 enddo ! End loop kdir 679 680! Increment by second vector 681 jj = jj + jdir 682 rxj = rxj + rcx2 683 ryj = ryj + rcy2 684 rzj = rzj + rcz2 685! Check to see if this direction is complete 686 lallfound2 = (r2.gt.r2j.and.r2.gt.MMcutoff2.and. & 687 nvec == nveck0) 688 r2j = r2 689 enddo ! End loop lallfound2 690 enddo ! End loop jdir 691 692! Increment by first vector 693 ii = ii + idir 694 rxi = rxi + rcx1 695 ryi = ryi + rcy1 696 rzi = rzi + rcz1 697 698! Check to see if this direction is complete 699 lallfound1 = (r2.gt.r2i.and.r2.gt.MMcutoff2 .and. & 700 nvec == nvecj0) 701 r2i = r2 702 enddo ! End loop lallfound1 703 enddo ! End loop idir 704 endif ! Endif over valid potentials 705 enddo ! End loops over atoms 706 enddo 707 708! 709! Print and add MM contribution to stress 710! 711 if (istr.ne.0) then 712 if (IONode .and. nMMpot > 0) then 713 write(6,'(/,a,6f12.2)') 'MM-Stress (kbar):', & 714 (mm_stress(jx,jx)/kbar,jx=1,3), & 715 mm_stress(1,2)/kbar, & 716 mm_stress(2,3)/kbar, & 717 mm_stress(1,3)/kbar 718 endif 719 720 stress = stress + mm_stress 721 endif 722 723! 724! Free workspace arrays 725! 726 call de_alloc( lvalidpot, 'lvalidpot', 'twobody' ) 727! 728! Stop timer 729! 730 call timer('MolMec', 2 ) 731 732 end subroutine twobody 733 734 735 subroutine plot_functions() 736! 737! Writes out V(r) info to files of the form MMpot.NN 738! Units: energy: eV, distance: Ang 739! 740 use units, only : eV, Ang 741 742 integer :: np 743 character(len=20) :: fname 744 integer, parameter :: npts = 1000 745 746 real(dp) :: rmin, rmax, delta, range 747 real(dp) :: etrm, ftrm, r1, r2, factor, br6, ebr6, f6, etrm1 748 real(dp) :: df6, ebr8, br8, f8, df8, br10, ebr10, f10, df10, r 749 real(dp) :: fg, fgprime, arg, earg, R0 750 integer :: i, iu 751 752 factor = 1.0_dp !! ?? 753 754 do np = 1,nMMpot 755 write(fname,"(a,i2.2)") "MMpot.", np 756 call io_assign(iu) 757 open(iu,file=trim(fname),form="formatted",status="replace", & 758 position="rewind",action="write") 759 760 rmin = 0.1_dp 761 rmax = min(20.0_dp,MMcutoff) 762 range = rmax - rmin 763 delta = range/npts 764 do i = 0, npts 765 r1 = rmin + delta * i 766 r2 = r1*r1 767 768 select case ( MMpottype(np) ) 769 case ( 1 ) ! C6 770 if (MMpotpar(2,np) == 0.0_dp) then 771 etrm = MMpotpar(1,np)/(r2**3) 772 ftrm = 6.0_dp*factor*etrm/r2 773 etrm = - etrm 774 else 775 br6 = MMpotpar(2,np)*sqrt(r2) 776 ebr6 = exp(-br6) 777 f6 = 1.0_dp + br6*(1.0_dp + 0.5_dp*br6*(1.0_dp + (br6/3.0_dp)*( & 778 1.0_dp + 0.25_dp*br6*(1.0_dp + 0.2_dp*br6*(1.0_dp + (br6/6.0_dp)))))) 779 f6 = 1.0_dp - f6*ebr6 780 etrm1 = MMpotpar(1,np)/(r2**3) 781 df6 = ebr6*(br6*(br6**6))/720.0_dp 782 etrm = - etrm1*f6 783 ftrm = factor*(6.0_dp*etrm1*f6 - etrm1*df6)/r2 784 endif 785 case ( 2 ) ! C8 786 if (MMpotpar(2,np) == 0.0_dp) then 787 etrm = MMpotpar(1,np)/(r2**4) 788 ftrm = 8.0_dp*factor*etrm/r2 789 etrm = - etrm 790 else 791 br8 = MMpotpar(2,np)*sqrt(r2) 792 ebr8 = exp(-br8) 793 f8 = 1.0_dp + br8*(1.0_dp + 0.5_dp*br8*(1.0_dp + (br8/3.0_dp)*( & 794 1.0_dp + 0.25_dp*br8*(1.0_dp + 0.2_dp*br8*(1.0_dp + & 795 (br8/6.0_dp)*(1.0_dp+(br8/7.0_dp)*(1.0_dp + 0.125_dp*br8))))))) 796 f8 = 1.0_dp - f8*ebr8 797 etrm1 = MMpotpar(1,np)/(r2**4) 798 df8 = ebr8*(br8*(br8**8))/40320.0_dp 799 etrm = - etrm1*f8 800 ftrm = factor*(8.0_dp*etrm1*f8 - etrm1*df8)/r2 801 endif 802 case ( 3 ) ! C10 803 if (MMpotpar(2,np) == 0.0_dp) then 804 etrm = MMpotpar(1,np)/(r2**5) 805 ftrm = 10.0_dp*factor*etrm/r2 806 etrm = - etrm 807 else 808 br10 = MMpotpar(2,np)*sqrt(r2) 809 ebr10 = exp(-br10) 810 f10 = 1.0_dp + br10*(1.0_dp + 0.5_dp*br10*(1.0_dp + & 811 (br10/3.0_dp)*(1.0_dp + 0.25_dp*br10*(1.0_dp + 0.2_dp*br10*( & 812 1.0_dp + (br10/6.0_dp)*(1.0_dp + (br10/7.0_dp)*(1.0_dp + & 813 0.125_dp*br10*(1.0_dp + (br10/9.0_dp)*(1.0_dp + 0.1_dp*br10))))))))) 814 f10 = 1.0_dp - f10*ebr10 815 etrm1 = MMpotpar(1,np)/(r2**5) 816 df10 = ebr10*(br10*(br10**10))/3628800.0_dp 817 etrm = - etrm1*f10 818 ftrm = factor*(10.0_dp*etrm1*f10 - etrm1*df10)/r2 819 endif 820 case ( 4 ) ! HARM 821 r = sqrt(r2) 822 etrm = MMpotpar(1,np)*(r - MMpotpar(2,np)) 823 ftrm = factor*etrm/r 824 etrm = 0.5_dp*etrm*(r - MMpotpar(2,np)) 825 case ( 5 ) ! Grimme 826 r = sqrt(r2) 827 R0 = MMpotpar(2,np) 828 arg = - d_grimme*(r/R0 - 1.0_dp) 829 earg = exp(arg) 830 fg = 1.0_dp / ( 1.0_dp + earg ) 831 etrm1 = s6_grimme * MMpotpar(1,np)/(r2**3) 832 fgprime = (d_grimme/R0) * earg / ( 1 + earg )**2 833 etrm = - etrm1*fg 834 ftrm = factor*(6.0_dp*etrm1*fg - etrm1*r*fgprime)/r2 835 end select 836 837 write(iu,*) r1/Ang, etrm/eV, ftrm*Ang/eV 838 839 enddo !! points 840 call io_close(iu) 841 enddo !! nPots 842 843 end subroutine plot_functions 844 845 end module molecularmechanics 846