1!--------------------------------------------------------------------------------------------------! 2! FES: a fast and general program to map metadynamics on grids ! 3! Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013 ! 4! Teodoro Laino ! 5!-----------------------------------------------------------------------------! 6 7! ************************************************************************************************** 8!> \brief Program to Map on grid the hills spawned during a metadynamics run 9!> \author Teodoro Laino [tlaino] - 06.2009 10!> \par History 11!> 03.2006 created [tlaino] 12!> teodoro.laino .at. gmail.com 13!> 11.2007 - tlaino (University of Zurich): Periodic COLVAR - cleaning. 14!> 12.2010 - teodoro.laino@gmail.com: addition of the MEP for FES 15!> 16!> \par Note 17!> Please report any bug to the author 18! ************************************************************************************************** 19MODULE graph_methods 20 21 USE cp_files, ONLY: close_file,& 22 open_file 23 USE graph_utils, ONLY: derivative,& 24 get_val_res,& 25 mep_input_data_type,& 26 pbc,& 27 point_no_pbc,& 28 point_pbc 29 USE kinds, ONLY: dp 30 USE memory_utilities, ONLY: reallocate 31 USE periodic_table, ONLY: init_periodic_table,& 32 nelem,& 33 ptable 34 USE physcon, ONLY: bohr 35 USE string_utilities, ONLY: uppercase 36#include "../base/base_uses.f90" 37 38 IMPLICIT NONE 39 PRIVATE 40 41 PUBLIC :: fes_compute_low, & 42 fes_write, & 43 fes_only_write, & 44 fes_min, & 45 fes_path, & 46 fes_cube_write 47 48CONTAINS 49! ************************************************************************************************** 50!> \brief Efficiently map the gaussians on the grid 51!> \param idim ... 52!> \param nn ... 53!> \param fes ... 54!> \param gauss ... 55!> \param ind ... 56!> \param ind0 ... 57!> \param nfes ... 58!> \param ndim ... 59!> \param ngauss ... 60!> \param ngrid ... 61!> \param iperd ... 62!> \par History 63!> 03.2006 created [tlaino] 64!> teodoro.laino .at. gmail.com 65!> \author Teodoro Laino 66! ************************************************************************************************** 67 RECURSIVE SUBROUTINE fes_compute_low(idim, nn, fes, gauss, ind, ind0, nfes, ndim, & 68 ngauss, ngrid, iperd) 69 INTEGER, INTENT(in) :: idim 70 INTEGER, DIMENSION(:) :: nn 71 REAL(KIND=dp), DIMENSION(:), POINTER :: fes 72 REAL(KIND=dp), DIMENSION(:, :), POINTER :: gauss 73 INTEGER, DIMENSION(:) :: ind, ind0 74 INTEGER, INTENT(in) :: nfes, ndim, ngauss 75 INTEGER, DIMENSION(:), POINTER :: ngrid 76 INTEGER, DIMENSION(:) :: iperd 77 78 INTEGER :: i, j, k, pnt 79 INTEGER, DIMENSION(:), POINTER :: ll, pos 80 REAL(KIND=dp) :: prod 81 82 ALLOCATE (pos(ndim), ll(ndim)) 83 pos = ind 84 k = nn(idim) 85 86 DO i = -k, k 87 pos(idim) = ind(idim) + i 88 IF (iperd(idim) == 0) THEN 89 IF (pos(idim) .GT. ngrid(idim)) CYCLE 90 IF (pos(idim) .LT. 1) CYCLE 91 END IF 92 IF (idim /= 1) THEN 93 CALL fes_compute_low(idim - 1, nn, fes, gauss, pos, ind0, nfes, ndim, ngauss, ngrid, iperd) 94 ELSE 95 pnt = point_pbc(pos, iperd, ngrid, ndim) 96 prod = 1.0_dp 97 DO j = 1, ndim 98 ll(j) = pos(j) - ind0(j) 99 prod = prod*gauss(ll(j), j) 100 END DO 101 fes(pnt) = fes(pnt) + prod 102 END IF 103 END DO 104 DEALLOCATE (pos, ll) 105 106 END SUBROUTINE fes_compute_low 107 108! ************************************************************************************************** 109!> \brief Writes the FES on the file 110!> \param unit_nr ... 111!> \param idim ... 112!> \param fes ... 113!> \param pos ... 114!> \param ndim ... 115!> \param ngrid ... 116!> \param dp_grid ... 117!> \param x0 ... 118!> \param ndw ... 119!> \param l_fes_int ... 120!> \param array ... 121!> \par History 122!> 03.2006 created [tlaino] 123!> teodoro.laino .at. gmail.com 124!> \author Teodoro Laino 125! ************************************************************************************************** 126 RECURSIVE SUBROUTINE fes_write(unit_nr, idim, fes, pos, ndim, ngrid, & 127 dp_grid, x0, ndw, l_fes_int, array) 128 INTEGER, INTENT(IN) :: unit_nr, idim 129 REAL(KIND=dp), DIMENSION(:), POINTER :: fes 130 INTEGER, DIMENSION(:), POINTER :: pos 131 INTEGER, INTENT(IN) :: ndim 132 INTEGER, DIMENSION(:), POINTER :: ngrid 133 REAL(KIND=dp), DIMENSION(:), POINTER :: dp_grid, x0 134 INTEGER, INTENT(IN) :: ndw 135 LOGICAL, INTENT(IN) :: l_fes_int 136 REAL(KIND=dp), DIMENSION(:), OPTIONAL :: array 137 138 INTEGER :: dimval, i, id, ind, is, it, itt, np, pnt 139 REAL(KIND=dp) :: dvol, sum_fes 140 REAL(KIND=dp), DIMENSION(:), POINTER :: xx 141 142 ALLOCATE (xx(ndim)) 143 xx = x0 144 DO i = 1, ngrid(idim) 145 pos(idim) = i 146 IF (idim /= ndim - ndw + 1) THEN 147 IF (PRESENT(array)) THEN 148 CALL fes_write(unit_nr, idim - 1, fes, pos, ndim, ngrid, dp_grid, & 149 x0, ndw, l_fes_int, array) 150 ELSE 151 CALL fes_write(unit_nr, idim - 1, fes, pos, ndim, ngrid, dp_grid, & 152 x0, ndw, l_fes_int) 153 END IF 154 ELSE 155 IF (PRESENT(array)) THEN 156 ind = 1 157 np = ngrid(ndim)*ngrid(ndim - 1)*ngrid(ndim - 2) 158 DO is = 1, ndw 159 itt = 1 160 DO it = 1, is - 1 161 itt = itt*ngrid(ndim - it) 162 END DO 163 ind = ind + (pos(ndim - is + 1) - 1)*itt 164 END DO 165 IF (ind > np) CPABORT("something wrong in indexing ..") 166 END IF 167 pnt = point_no_pbc(pos, ngrid, ndim) 168 xx = x0 + dp_grid*(pos - 1) 169 dimval = PRODUCT(ngrid(1:ndim - ndw)) 170 171 IF (.NOT. l_fes_int) THEN 172 IF (PRESENT(array)) THEN 173 array(ind) = MINVAL(-fes(pnt:pnt + dimval - 1)) 174 ELSE 175 WRITE (unit_nr, '(10f20.10)') (xx(id), id=ndim, ndim - ndw + 1, -1), MINVAL(-fes(pnt:pnt + dimval - 1)) 176 END IF 177 ELSE 178 sum_fes = 0.0_dp 179 dvol = 1.0_dp 180 dvol = PRODUCT(dp_grid(1:ndim - ndw)) 181 DO is = pnt, pnt + dimval - 1 182 sum_fes = sum_fes + fes(is)*dvol 183 END DO 184 IF (PRESENT(array)) THEN 185 array(ind) = -sum_fes 186 ELSE 187 WRITE (unit_nr, '(10f20.10)') (xx(id), id=ndim, ndim - ndw + 1, -1), -sum_fes 188 END IF 189 END IF 190 END IF 191 END DO 192 DEALLOCATE (xx) 193 194 END SUBROUTINE fes_write 195 196! ************************************************************************************************** 197!> \brief Writes the FES on the file when stride is requested 198!> \param idim ... 199!> \param fes ... 200!> \param pos ... 201!> \param ndim ... 202!> \param ngrid ... 203!> \param dp_grid ... 204!> \param ndw ... 205!> \param l_fes_int ... 206!> \param unit_nr ... 207!> \par History 208!> 03.2006 created [tlaino] 209!> teodoro.laino .at. gmail.com 210!> \author Teodoro Laino 211! ************************************************************************************************** 212 RECURSIVE SUBROUTINE fes_only_write(idim, fes, pos, ndim, ngrid, dp_grid, ndw, l_fes_int, unit_nr) 213 INTEGER, INTENT(IN) :: idim 214 REAL(KIND=dp), DIMENSION(:), POINTER :: fes 215 INTEGER, DIMENSION(:), POINTER :: pos 216 INTEGER, INTENT(IN) :: ndim 217 INTEGER, DIMENSION(:), POINTER :: ngrid 218 REAL(KIND=dp), DIMENSION(:), POINTER :: dp_grid 219 INTEGER, INTENT(IN) :: ndw 220 LOGICAL, INTENT(IN) :: l_fes_int 221 INTEGER :: unit_nr 222 223 INTEGER :: dimval, i, is, pnt 224 REAL(KIND=dp) :: dvol, sum_fes 225 226 DO i = 1, ngrid(idim) 227 pos(idim) = i 228 IF (idim /= ndim - ndw + 1) THEN 229 CALL fes_only_write(idim - 1, fes, pos, ndim, ngrid, dp_grid, ndw, l_fes_int, unit_nr) 230 ELSE 231 pnt = point_no_pbc(pos, ngrid, ndim) 232 dimval = PRODUCT(ngrid(1:ndim - ndw)) 233 IF (l_fes_int) THEN 234 WRITE (unit_nr, '(1f12.5)') MINVAL(-fes(pnt:pnt + dimval - 1)) 235 ELSE 236 sum_fes = 0.0_dp 237 dvol = PRODUCT(dp_grid(1:ndim - ndw)) 238 DO is = pnt, pnt + dimval - 1 239 sum_fes = sum_fes + fes(is)*dvol 240 END DO 241 WRITE (unit_nr, '(1f12.5)') - sum_fes 242 END IF 243 END IF 244 END DO 245 246 END SUBROUTINE fes_only_write 247 248! ************************************************************************************************** 249!> \brief Finds minima of the FES 250!> \param fes ... 251!> \param ndim ... 252!> \param iperd ... 253!> \param ngrid ... 254!> \param dp_grid ... 255!> \param x0 ... 256!> \param ndw ... 257!> \par History 258!> 06.2009 created [tlaino] 259!> teodoro.laino .at. gmail.com 260!> \author Teodoro Laino 261! ************************************************************************************************** 262 SUBROUTINE fes_min(fes, ndim, iperd, ngrid, dp_grid, x0, ndw) 263 REAL(KIND=dp), DIMENSION(:), POINTER :: fes 264 INTEGER, INTENT(IN) :: ndim 265 INTEGER, DIMENSION(:), POINTER :: iperd, ngrid 266 REAL(KIND=dp), DIMENSION(:), POINTER :: dp_grid, x0 267 INTEGER, INTENT(IN) :: ndw 268 269 INTEGER :: i, id, iter, j, k, max_ntrust, nacc, & 270 ntrials, pnt 271 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: history 272 INTEGER, DIMENSION(:), POINTER :: pos, pos0 273 INTEGER, DIMENSION(ndim) :: Dpos, ntrust 274 LOGICAL :: do_save 275 REAL(KIND=dp) :: fes_now, fes_old, norm_dx, resto 276 REAL(KIND=dp), DIMENSION(:), POINTER :: dx, rnd, xx 277 278 IF (ndw /= ndim) CPABORT("Not implemented for projected FES!") 279 280 ntrust = ngrid/10 281 ntrials = PRODUCT(ngrid) 282 WRITE (*, '(A,10I6)', ADVANCE="no") "FES| Trust hyper-radius ", ntrust 283 WRITE (*, '(A,10F12.6)') " which is equivalent to: ", ntrust*dp_grid 284 285 ALLOCATE (xx(ndim), dx(ndim), pos0(ndim), rnd(ndim), pos(ndim)) 286 ALLOCATE (history(ndim, ntrials)) 287 history = 0 288 nacc = 0 289 Trials: DO j = 1, ntrials 290 ! Loop over all points 291 pnt = j 292 DO k = ndim, 2, -1 293 pos0(k) = pnt/PRODUCT(ngrid(1:k - 1)) 294 resto = MOD(pnt, PRODUCT(ngrid(1:k - 1))) 295 IF (resto /= 0) THEN 296 pnt = pnt - pos0(k)*PRODUCT(ngrid(1:k - 1)) 297 pos0(k) = pos0(k) + 1 298 ELSE 299 pnt = PRODUCT(ngrid(1:k - 1)) 300 END IF 301 END DO 302 pos0(1) = pnt 303 304 ! Loop over the frame points unless it is periodic 305 DO k = 1, ndim 306 IF ((iperd(k) == 0) .AND. (pos0(k) < ntrust(k))) CYCLE Trials 307 IF ((iperd(k) == 0) .AND. (pos0(k) > ngrid(k) - ntrust(k))) CYCLE Trials 308 END DO 309 310 ! Evaluate position and derivative 311 pos = pos0 312 xx = x0 + dp_grid*(pos - 1) 313 dx = derivative(fes, pos, iperd, ndim, ngrid, dp_grid) 314 315 ! Integrate till derivative is small enough.. 316 pnt = point_no_pbc(pos, ngrid, ndim) 317 fes_now = -fes(pnt) 318 fes_old = HUGE(0.0_dp) 319 320 i = 1 321 DO WHILE ((i <= 100) .OR. (fes_now < fes_old)) 322 fes_old = fes_now 323 !WRITE(10+j,'(10f20.10)')(xx(id),id=ndim,1,-1),-fes(pnt) 324 325 norm_dx = SQRT(DOT_PRODUCT(dx, dx)) 326 IF (norm_dx == 0.0_dp) EXIT ! It is in a really flat region 327 xx = xx - MIN(0.1_dp, norm_dx)*dx/norm_dx 328 ! Re-evaluating pos 329 pos = CEILING((xx - x0)/dp_grid) + 1 330 CALL pbc(pos, iperd, ngrid, ndim) 331 332 ! Incremental pos 333 dx = derivative(fes, pos, iperd, ndim, ngrid, dp_grid) 334 pnt = point_no_pbc(pos, ngrid, ndim) 335 fes_now = -fes(pnt) 336 i = i + 1 337 END DO 338 iter = i 339 340 ! Compare with the available minima and if they are the same skip 341 ! saving this position.. 342 do_save = fes(pnt) >= 1.0E-3_dp 343 IF (do_save) THEN 344 DO i = 1, nacc 345 Dpos = pos - history(:, i) 346 norm_dx = DOT_PRODUCT(Dpos, Dpos) 347 max_ntrust = MAXVAL(ntrust) 348 ! (SQRT(REAL(norm_dx, KIND=dp)) <= MAXVAL(ntrust)) ... 349 IF ((norm_dx <= REAL(max_ntrust*max_ntrust, KIND=dp)) .OR. (fes(pnt) < 1.0E-3_dp)) THEN 350 do_save = .FALSE. 351 EXIT 352 END IF 353 END DO 354 END IF 355 IF (do_save) THEN 356 pnt = point_no_pbc(pos, ngrid, ndim) 357 xx = x0 + dp_grid*(pos - 1) 358 WRITE (*, '(A,5F12.6)', ADVANCE="NO") "FES| Minimum found (", (xx(id), id=ndim, ndim - ndw + 1, -1) 359 WRITE (*, '(A,F12.6,A,I6)') " ). FES value = ", -fes(pnt), " Hartree. Number of Iter: ", iter 360 nacc = nacc + 1 361 history(:, nacc) = pos 362 END IF 363 END DO Trials 364 WRITE (*, '(A,I6,A)') "FES| Number of Minimum found: ", nacc, "." 365 366 DEALLOCATE (xx, dx, pos0, rnd, pos) 367 DEALLOCATE (history) 368 369 END SUBROUTINE fes_min 370 371! ************************************************************************************************** 372!> \brief Finds path between two points (a) and (b) 373!> \param fes ... 374!> \param ndim ... 375!> \param ngrid ... 376!> \param dp_grid ... 377!> \param iperd ... 378!> \param x0 ... 379!> \param ndw ... 380!> \param mep_input_data ... 381!> \param l_int ... 382!> \par History 383!> 12.2010 created [tlaino] 384!> teodoro.laino .at. gmail.com 385!> \author Teodoro Laino 386! ************************************************************************************************** 387 SUBROUTINE fes_path(fes, ndim, ngrid, dp_grid, iperd, x0, ndw, mep_input_data, l_int) 388 REAL(KIND=dp), DIMENSION(:), POINTER :: fes 389 INTEGER, INTENT(IN) :: ndim 390 INTEGER, DIMENSION(:), POINTER :: ngrid 391 REAL(KIND=dp), DIMENSION(:), POINTER :: dp_grid 392 INTEGER, DIMENSION(:), POINTER :: iperd 393 REAL(KIND=dp), DIMENSION(:), POINTER :: x0 394 INTEGER, INTENT(IN) :: ndw 395 TYPE(mep_input_data_type), INTENT(IN) :: mep_input_data 396 LOGICAL :: l_int 397 398 INTEGER :: i, id, irep, iter, nf, nreplica, ns, & 399 pnt, unit_nr 400 INTEGER, DIMENSION(:), POINTER :: ipos 401 LOGICAL :: converged 402 REAL(KIND=dp) :: avg1, avg2, diff, ene, norm_dx, xx0, yy0 403 REAL(KIND=dp), DIMENSION(:), POINTER :: davg1, davg2, dxx, dyy, fes_rep, tang, & 404 xx, yy 405 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dx, pos, pos_old 406 407 IF (ndw /= ndim) CPABORT("Not implemented for projected FES!") 408 nreplica = mep_input_data%nreplica 409 ALLOCATE (xx(ndim), dx(ndim, nreplica), pos_old(ndim, nreplica), pos(ndim, nreplica), & 410 ipos(ndim), fes_rep(nreplica), dxx(ndim), dyy(ndim), yy(ndim), davg1(ndim), & 411 tang(ndim), davg2(ndim)) 412 413 IF (l_int) THEN 414 id = 0 415 DO i = ndim, ndim - ndw + 1, -1 416 id = id + 1 417 pos(i, 1) = mep_input_data%minima(id, 1) 418 pos(i, nreplica) = mep_input_data%minima(id, 2) 419 END DO 420 421 ! Interpolate nreplica-2 points 422 xx = (pos(:, nreplica) - pos(:, 1))/REAL(nreplica - 1, KIND=dp) 423 DO irep = 2, nreplica - 1 424 pos(:, irep) = pos(:, 1) + xx(:)*REAL(irep - 1, KIND=dp) 425 END DO 426 427 ELSE 428 pos = mep_input_data%minima 429 ENDIF 430 431 ! Compute value and derivative in all replicas 432 DO irep = 1, nreplica 433 ipos = FLOOR((pos(:, irep) - x0)/dp_grid) + 1 434 pnt = point_no_pbc(ipos, ngrid, ndim) 435 dx(:, irep) = derivative(fes, ipos, iperd, ndim, ngrid, dp_grid) 436 fes_rep(irep) = -fes(pnt) 437 END DO 438 439 ! Implement a simple elastic band method (Hamiltonian): definitely not the best 440 ! method, but for such a simple task it should be more than enough 441 converged = .FALSE. 442 pos_old = pos 443 iter = 0 444 DO WHILE ((.NOT. converged) .AND. (iter <= mep_input_data%max_iter)) 445 iter = iter + 1 446 avg1 = 0.0_dp 447 ! compute average length (distance 1) 448 DO irep = 2, nreplica 449 xx = pos(:, irep) - pos(:, irep - 1) 450 avg1 = avg1 + SQRT(DOT_PRODUCT(xx, xx)) 451 END DO 452 avg1 = avg1/REAL(nreplica - 1, KIND=dp) 453 454 avg2 = 0.0_dp 455 ! compute average length (distance 2) 456 DO irep = 3, nreplica 457 xx = pos(:, irep) - pos(:, irep - 2) 458 avg2 = avg2 + SQRT(DOT_PRODUCT(xx, xx)) 459 END DO 460 avg2 = avg2/REAL(nreplica - 2, KIND=dp) 461 462 ! compute energy and derivatives 463 dx = 0.0_dp 464 ene = 0.0_dp 465 ns = 1 466 nf = nreplica 467 DO irep = 1, nreplica 468 ! compute energy and map point replica irep 469 ipos = FLOOR((pos(:, irep) - x0)/dp_grid) + 1 470 pnt = point_no_pbc(ipos, ngrid, ndim) 471 fes_rep(irep) = -fes(pnt) 472 IF ((irep == 1) .OR. (irep == nreplica)) CYCLE 473 474 ! ------------------------------------------------------------- 475 ! compute non-linear elastic terms : including only 2-d springs 476 ! ------------------------------------------------------------- 477 davg2 = 0.0_dp 478 IF (irep < nf - 1) THEN 479 xx = pos(:, irep) - pos(:, irep + 2) 480 xx0 = SQRT(DOT_PRODUCT(xx, xx)) 481 dxx = 1.0_dp/xx0*xx 482 ene = ene + 0.25_dp*mep_input_data%kb*(xx0 - avg2)**2 483 davg2 = davg2 + dxx 484 END IF 485 486 IF (irep > ns + 1) THEN 487 xx = pos(:, irep) - pos(:, irep - 2) 488 yy0 = SQRT(DOT_PRODUCT(xx, xx)) 489 dyy = 1.0_dp/yy0*xx 490 davg2 = davg2 + dyy 491 END IF 492 davg2 = davg2/REAL(nreplica - 2, KIND=dp) 493 494 IF (irep < nf - 1) THEN 495 dx(:, irep) = dx(:, irep) + 0.5_dp*mep_input_data%kb*(xx0 - avg2)*(dxx - davg2) 496 END IF 497 IF (irep > ns + 1) THEN 498 dx(:, irep) = dx(:, irep) + 0.5_dp*mep_input_data%kb*(yy0 - avg2)*(dyy - davg2) 499 END IF 500 501 ! ------------------------------------------------------------- 502 ! Evaluation of the elastic term 503 ! ------------------------------------------------------------- 504 xx = pos(:, irep) - pos(:, irep + 1) 505 yy0 = SQRT(DOT_PRODUCT(xx, xx)) 506 dyy = 1.0_dp/yy0*xx 507 508 xx = pos(:, irep) - pos(:, irep - 1) 509 xx0 = SQRT(DOT_PRODUCT(xx, xx)) 510 dxx = 1.0_dp/xx0*xx 511 davg1 = (dxx + dyy)/REAL(nreplica - 1, KIND=dp) 512 513 ene = ene + 0.5_dp*mep_input_data%kb*(xx0 - avg1)**2 514 dx(:, irep) = dx(:, irep) + mep_input_data%kb*(xx0 - avg1)*(dxx - davg1) + & 515 mep_input_data%kb*(yy0 - avg1)*(dyy - davg1) 516 517 ! Evaluate the tangent 518 xx = pos(:, irep + 1) - pos(:, irep) 519 xx = xx/SQRT(DOT_PRODUCT(xx, xx)) 520 yy = pos(:, irep) - pos(:, irep - 1) 521 yy = yy/SQRT(DOT_PRODUCT(yy, yy)) 522 tang = xx + yy 523 tang = tang/SQRT(DOT_PRODUCT(tang, tang)) 524 525 xx = derivative(fes, ipos, iperd, ndim, ngrid, dp_grid) 526 dx(:, irep) = DOT_PRODUCT(dx(:, irep), tang)*tang + & 527 xx - DOT_PRODUCT(xx, tang)*tang 528 END DO 529 dx(:, 1) = 0.0_dp 530 dx(:, nreplica) = 0.0_dp 531 532 ! propagate the band with a SD step 533 diff = 0.0_dp 534 DO irep = 1, nreplica 535 ene = ene + fes_rep(irep) 536 IF ((irep == 1) .OR. (irep == nreplica)) CYCLE 537 538 norm_dx = SQRT(DOT_PRODUCT(dx(:, irep), dx(:, irep))) 539 IF (norm_dx /= 0.0_dp) THEN 540 pos(:, irep) = pos(:, irep) - MIN(0.1_dp, norm_dx)*dx(:, irep)/norm_dx 541 END IF 542 xx = pos(:, irep) - pos_old(:, irep) 543 diff = diff + DOT_PRODUCT(xx, xx) 544 END DO 545 ! SQRT(diff) <= 0.001_dp 546 IF (diff <= 1.0e-6_dp) THEN 547 converged = .TRUE. 548 END IF 549 pos_old = pos 550 WRITE (*, *) "Iteration nr.", iter, SQRT(diff) 551 END DO 552 553 WRITE (*, *) "MEP saved on <mep.data> file." 554 CALL open_file(unit_number=unit_nr, file_name="mep.data", file_action="WRITE", file_status="UNKNOWN", file_form="FORMATTED") 555 DO irep = 1, nreplica 556 ! compute energy and derivative for each single point of the replica 557 ipos = FLOOR((pos(:, irep) - x0)/dp_grid) + 1 558 pnt = point_no_pbc(ipos, ngrid, ndim) 559 fes_rep(irep) = -fes(pnt) 560 WRITE (unit_nr, *) irep, pos(:, nreplica - irep + 1), fes_rep(nreplica - irep + 1) 561 END DO 562 CALL close_file(unit_nr) 563 564 DEALLOCATE (xx, dx, pos, fes_rep, ipos, pos_old, yy, dyy, dxx, davg1, tang, davg2) 565 END SUBROUTINE fes_path 566 567! ************************************************************************************************** 568!> \brief Dump FES with a GAUSSIAN cube format - Useful for multidimensional FES 569!> \param idim ... 570!> \param fes ... 571!> \param pos ... 572!> \param ndim ... 573!> \param ngrid ... 574!> \param dp_grid ... 575!> \param x0 ... 576!> \param ndw ... 577!> \param l_fes_int ... 578!> \param file ... 579!> \par History 580!> 12.2013 created [tlaino] 581!> teodoro.laino .at. gmail.com 582!> \author Teodoro Laino 583! ************************************************************************************************** 584 SUBROUTINE fes_cube_write(idim, fes, pos, ndim, ngrid, dp_grid, x0, ndw, l_fes_int, file) 585 INTEGER, INTENT(IN) :: idim 586 REAL(KIND=dp), DIMENSION(:), POINTER :: fes 587 INTEGER, DIMENSION(:), POINTER :: pos 588 INTEGER, INTENT(IN) :: ndim 589 INTEGER, DIMENSION(:), POINTER :: ngrid 590 REAL(KIND=dp), DIMENSION(:), POINTER :: dp_grid, x0 591 INTEGER, INTENT(IN) :: ndw 592 LOGICAL, INTENT(IN) :: l_fes_int 593 CHARACTER(LEN=80) :: file 594 595 CHARACTER(LEN=120) :: line 596 CHARACTER(LEN=5) :: label, labelp 597 INTEGER :: i, id(3), ii, iix, iiy, iiz, ix, iy, iz, & 598 natoms, np 599 INTEGER, DIMENSION(:), POINTER :: izat 600 REAL(KIND=dp) :: cell(3, 3), delta(3), dr(3), residual, & 601 rt(3) 602 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: rho, rhot 603 REAL(KIND=dp), DIMENSION(:, :), POINTER :: xat 604 605 CALL init_periodic_table() 606 IF (ndw .GT. 3) THEN 607 WRITE (*, *) 608 WRITE (*, *) 'ERROR: GAUSSIAN format can only handle FES on 3 CV !' 609 CPABORT("") 610 END IF 611 612 OPEN (10, file=file, status='old') 613 CALL get_val_res(unit=10, section="&SUBSYS", subsection="&CELL") 614 READ (10, *) label, cell(1, 1), cell(2, 1), cell(3, 1) 615 READ (10, *) label, cell(1, 2), cell(2, 2), cell(3, 2) 616 READ (10, *) label, cell(1, 3), cell(2, 3), cell(3, 3) 617 rt(1) = -(cell(1, 1)/2._dp) 618 rt(2) = -(cell(2, 2)/2._dp) 619 rt(3) = -(cell(3, 3)/2._dp) 620 621 WRITE (*, *) 'Dumping GAUSSIAN CUBE format' 622 WRITE (*, *) 'Cell vectors' 623 WRITE (*, *) 624 625 residual = 0.0d0 626 DO ix = 1, 3 627 DO iy = ix + 1, 3 628 residual = residual + cell(ix, iy)**2 629 END DO 630 END DO 631 632 IF (residual .GT. 1.0d-6) THEN 633 WRITE (*, *) 634 WRITE (*, *) 'ERROR: this program can only handle orthogonal cells' 635 WRITE (*, *) ' with vectors pointing in the X, Y and Z directions' 636 CPABORT("") 637 END IF 638 639 WRITE (*, *) 640 WRITE (*, *) 'Cube grid mesh: ', ngrid(1), 'x', ngrid(2), 'x', ngrid(3) 641 WRITE (*, *) 'Origin in:', rt 642 WRITE (*, *) 643 644 DO ix = 1, 3 645 dr(ix) = cell(ix, ix)/REAL(ngrid(ix), KIND=dp) 646 END DO 647 648 np = ngrid(1)*ngrid(2)*ngrid(3) 649 ALLOCATE (rho(np), rhot(np)) 650 CALL fes_write(123, idim, fes, pos, ndim, ngrid, dp_grid, x0, ndw, l_fes_int, rho) 651 WRITE (*, *) 'Internal FES transfer completed!' 652 653 ! translate cell 654 DO ix = 1, 3 655 delta(ix) = rt(ix)/dr(ix) 656 id(ix) = INT(delta(ix)) 657 delta(ix) = rt(ix) - id(ix)*dr(ix) 658 END DO 659 660 DO iz = 1, ngrid(3) 661 DO iy = 1, ngrid(2) 662 DO ix = 1, ngrid(1) 663 iix = ix + id(1) 664 iiy = iy + id(2) 665 iiz = iz + id(3) 666 IF (iix .LT. 1) iix = iix + ngrid(1) 667 IF (iiy .LT. 1) iiy = iiy + ngrid(2) 668 IF (iiz .LT. 1) iiz = iiz + ngrid(3) 669 IF (iix .GT. ngrid(1)) iix = iix - ngrid(1) 670 IF (iiy .GT. ngrid(2)) iiy = iiy - ngrid(2) 671 IF (iiz .GT. ngrid(3)) iiz = iiz - ngrid(3) 672 673 IF (iix .LT. 1) CPABORT("ix < 0") 674 IF (iiy .LT. 1) CPABORT("iy < 0") 675 IF (iiz .LT. 1) CPABORT("iz < 0") 676 IF (iix .GT. ngrid(1)) CPABORT("ix > cell") 677 IF (iiy .GT. ngrid(2)) CPABORT("iy > cell") 678 IF (iiz .GT. ngrid(3)) CPABORT("iz > cell") 679 i = ix + (iy - 1)*ngrid(1) + (iz - 1)*ngrid(1)*ngrid(2) 680 ii = iix + (iiy - 1)*ngrid(1) + (iiz - 1)*ngrid(1)*ngrid(2) 681 rhot(ii) = rho(i) 682 END DO 683 END DO 684 END DO 685 686 REWIND (10) 687 CALL get_val_res(unit=10, section="&SUBSYS", subsection="&COORD") 688 natoms = 0 689 ALLOCATE (xat(1000, 3)) 690 ALLOCATE (izat(1000)) 691 DO WHILE (.TRUE.) 692 READ (10, '(A)') line 693 IF (INDEX(line, '&END') /= 0) EXIT 694 natoms = natoms + 1 695 READ (line, *) label, (xat(natoms, ix), ix=1, 3) 696 IF (natoms == SIZE(xat, 1)) THEN 697 CALL reallocate(xat, 1, SIZE(xat, 1)*2, 1, 3) 698 CALL reallocate(izat, 1, SIZE(izat)*2) 699 END IF 700 CALL uppercase(label) 701 DO i = 1, nelem 702 labelp = ptable(i)%symbol 703 CALL uppercase(labelp) 704 IF (TRIM(label) == TRIM(labelp)) EXIT 705 END DO 706 IF (i == nelem + 1) THEN 707 WRITE (*, *) TRIM(label), "In line: ", line 708 CPABORT("Element not recognized!") 709 END IF 710 izat(natoms) = i 711 END DO 712 CALL reallocate(xat, 1, natoms, 1, 3) 713 CALL reallocate(izat, 1, natoms) 714 715 DO i = 1, natoms 716 DO ix = 1, 3 717 xat(i, ix) = xat(i, ix) + rt(ix) - delta(ix) 718 IF (xat(i, ix) .LT. rt(ix)) xat(i, ix) = xat(i, ix) + cell(ix, ix) 719 IF (xat(i, ix) .GT. -rt(ix)) xat(i, ix) = xat(i, ix) - cell(ix, ix) 720 END DO 721 END DO 722 723 WRITE (123, *) "FES on CUBE" 724 WRITE (123, *) "created by fes in CP2K" 725 WRITE (123, '(i5,3f12.6)') natoms, rt(1:3)*bohr 726 727 DO ix = 1, 3 728 ii = ngrid(ix) 729 WRITE (123, '(i5,4f12.6)') ii, (cell(ix, iy)/ii*bohr, iy=1, 3) 730 END DO 731 732 DO i = 1, natoms 733 WRITE (123, '(i5,4f12.6)') izat(i), 0.0, (xat(i, ix)*bohr, ix=1, 3) 734 END DO 735 736 DO ix = 1, ngrid(1) 737 DO iy = 1, ngrid(2) 738 739 WRITE (123, '(6e13.5)') (rhot(ix + (iy - 1)*ngrid(1) + (iz - 1)*ngrid(1)& 740 &*ngrid(2)), iz=1, ngrid(3)) 741 END DO 742 END DO 743 DEALLOCATE (xat, rho, rhot) 744 745 END SUBROUTINE fes_cube_write 746 747END MODULE graph_methods 748