1c $Id$ 2 3c ********************************************** 4c * * 5c * tamd_initialize * 6c * * 7c ********************************************** 8 subroutine tamd_initialize(rtdb) 9 implicit none 10 integer rtdb 11 12#include "bafdecls.fh" 13#include "btdb.fh" 14#include "errquit.fh" 15#include "stdio.fh" 16#include "tamd.fh" 17 18* **** local variables **** 19 integer taskid,MASTER 20 parameter (MASTER=0) 21 logical value,cnmeta_read,ok 22 integer d,i,j,n1,n2,s1,s2,yshift,ntype,seed 23 real*8 x,y,z,r,rmax,dr,pi,da,maxcoord,a,b,n,m,sigma,w,r0,gamma,k 24 real*8 boundary(2),temp,fac,sn2(12),mass,ztemp,zinitial(2),n12 25 character*80 rtdb_name,potential_filename 26 character*60 gauss_filename 27 character*4 celement1,celement2 28 character*500 eqnstring 29 30 31* **** external functions **** 32 logical control_Nose 33 real*8 lattice_unita,ion_Temperature,control_Nose_Tr 34 real*8 control_time_step 35 character*7 c_index_name 36 integer control_it_out 37 external control_Nose 38 external lattice_unita,ion_Temperature,control_Nose_Tr 39 external control_time_step 40 external c_index_name 41 external control_it_out 42 43* initialize variables to silence compilers 44 sigma = 0.0 45 w = 0.0 46 47 call Parallel_taskid(taskid) 48 metaprintcount = 0 49 rtdb_name = 'tamd_print' 50 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,maxmetaprintcount)) 51 > maxmetaprintcount=100 52 53 rtdb_name = 'tamd_nmeta' 54 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,nmeta)) nmeta = 0 55 56 tamdfound = (nmeta.gt.0) 57 58 if (tamdfound) then 59 60 61 rtdb_name = 'tamd_update' 62 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,maxmetacount)) 63 > maxmetacount = 1 64 65 rtdb_name = 'tamd_zfixed' 66 if (.not.btdb_get(rtdb,rtdb_name,mt_log,1,zfixed)) 67 > zfixed = .false. 68 69 rtdb_name = 'tamd_metacount' 70 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,metacount)) 71 > metacount = 0 72 73 rtdb_name = 'tamd_print_shift' 74 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,metarayshift)) 75 > metarayshift = maxmetaprintcount 76 77 rtdb_name = 'tamd_metaraycount' 78 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,metaraycount)) 79 > metaraycount = 0 80 81 rtdb_name = 'tamd_tempered' 82 if (.not.btdb_get(rtdb,rtdb_name,mt_dbl,1,dTtempered)) 83 > dTtempered = -1.0d0 84 85 rtdb_name = 'tamd_boundary' 86 if (.not.btdb_get(rtdb,rtdb_name,mt_dbl,2,boundary)) then 87 boundary(1) = 0.0d0 88 boundary(2) = 0.0d0 89 end if 90 91 rtdb_name = 'tamd_sn2-surface' 92 if (.not.btdb_get(rtdb,rtdb_name,mt_dbl,12,sn2)) then 93 call dcopy(12,0.0d0,0,sn2,1) 94 end if 95 96 value = .true. 97 if (nmeta.gt.0) then 98 rtdb_name = 'tamd_nindxmeta' 99 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,nindxmeta)) 100 > nindxmeta = 0 101 rtdb_name = 'tamd_nparammeta' 102 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,nparammeta)) 103 > nparammeta = 0 104 rtdb_name = 'tamd_nxmeta_all' 105 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,nxmeta_all)) 106 > nxmeta_all = 0 107 108* **** allocate and read nxmeta,sindxmeta,sparameta,and indxmeta **** 109 value = value.and. 110 > BA_alloc_get(mt_int,nmeta,'nxmeta', 111 > nxmeta(2),nxmeta(1)) 112 value = value.and. 113 > BA_alloc_get(mt_int,nmeta,'sindxmeta', 114 > sindxmeta(2),sindxmeta(1)) 115 value = value.and. 116 > BA_alloc_get(mt_int,nmeta,'sparammeta', 117 > sparammeta(2),sparammeta(1)) 118 value = value.and. 119 > BA_alloc_get(mt_int,nindxmeta,'indxmeta', 120 > indxmeta(2),indxmeta(1)) 121 122 rtdb_name = 'tamd_nxmeta' 123 if (.not.btdb_get(rtdb,rtdb_name,mt_int,nmeta, 124 > int_mb(nxmeta(1)))) then 125 do d=1,nmeta 126 int_mb(nxmeta(1)+d-1) = 0 127 end do 128 end if 129 rtdb_name = 'tamd_sindxmeta' 130 if (.not.btdb_get(rtdb,rtdb_name,mt_int,nmeta, 131 > int_mb(sindxmeta(1)))) then 132 do d=1,nmeta 133 int_mb(sindxmeta(1)+d-1) = 0 134 end do 135 end if 136 rtdb_name = 'tamd_sparammeta' 137 if (.not.btdb_get(rtdb,rtdb_name,mt_int,nmeta, 138 > int_mb(sparammeta(1)))) then 139 do d=1,nmeta 140 int_mb(sindxmeta(1)+d-1) = 0 141 end do 142 end if 143 rtdb_name = 'tamd_indxmeta' 144 if (.not.btdb_get(rtdb,rtdb_name,mt_int,nindxmeta, 145 > int_mb(indxmeta(1)))) then 146 do d=1,nindxmeta 147 int_mb(indxmeta(1)+d-1) = 0 148 end do 149 end if 150 151 152* ***** allocate and read pmeta,ameta,bmeta,parammeta, and ymeta **** 153 value = value.and. 154 > BA_alloc_get(mt_int,nmeta,'pmeta',pmeta(2),pmeta(1)) 155 value = value.and. 156 > BA_alloc_get(mt_dbl,nmeta,'ameta',ameta(2),ameta(1)) 157 value = value.and. 158 > BA_alloc_get(mt_dbl,nmeta,'bmeta',bmeta(2),bmeta(1)) 159 value = value.and. 160 > BA_alloc_get(mt_dbl,nparammeta,'parammeta', 161 > parammeta(2),parammeta(1)) 162 value = value.and. 163 > BA_alloc_get(mt_dbl,nxmeta_all, 164 > 'ymeta',ymeta(2),ymeta(1)) 165 166 rtdb_name = 'tamd_pmeta' 167 if (.not.btdb_get(rtdb,rtdb_name,mt_int,nmeta, 168 > int_mb(pmeta(1)))) then 169 do d=1,nmeta 170 int_mb(pmeta(1)+d-1) = 0 171 end do 172 end if 173 rtdb_name = 'tamd_ameta' 174 if (.not.btdb_get(rtdb,rtdb_name,mt_dbl,nmeta, 175 > dbl_mb(ameta(1)))) then 176 do d=1,nmeta 177 dbl_mb(ameta(1)+d-1) = 0.0d0 178 end do 179 end if 180 rtdb_name = 'tamd_bmeta' 181 if (.not.btdb_get(rtdb,rtdb_name,mt_dbl,nmeta, 182 > dbl_mb(bmeta(1)))) then 183 do d=1,nmeta 184 dbl_mb(bmeta(1)+d-1) = 0.0d0 185 end do 186 end if 187 rtdb_name = 'tamd_parammeta' 188 if (.not.btdb_get(rtdb,rtdb_name,mt_dbl,nparammeta, 189 > dbl_mb(parammeta(1)))) then 190 do d=1,nparammeta 191 dbl_mb(parammeta(1)+d-1) = 0.0d0 192 end do 193 end if 194 195* **** set default limits **** 196 x = lattice_unita(1,1) 197 y = lattice_unita(2,1) 198 z = lattice_unita(3,1) 199 r = dsqrt(x*x+y*y+z*z) 200 rmax = r 201 x = lattice_unita(1,2) 202 y = lattice_unita(2,2) 203 z = lattice_unita(3,2) 204 r = dsqrt(x*x+y*y+z*z) 205 if (r.gt.rmax) rmax = r 206 x = lattice_unita(1,3) 207 y = lattice_unita(2,3) 208 z = lattice_unita(3,3) 209 r = dsqrt(x*x+y*y+z*z) 210 if (r.gt.rmax) rmax = r 211 do d=1,nmeta 212 s1 = int_mb(sindxmeta(1) +d-1) 213 if ((int_mb(indxmeta(1)+s1).eq.1).or. 214 > (int_mb(indxmeta(1)+s1).eq.5).or. 215 > (int_mb(indxmeta(1)+s1).eq.13)) then 216 if (int_mb(pmeta(1)+d-1).lt.0) 217 > int_mb(pmeta(1)+d-1) = 0 218 if (dbl_mb(ameta(1)+d-1).lt.0.0d0) 219 > dbl_mb(ameta(1)+d-1) = 0.0d0 220 if (dbl_mb(bmeta(1)+d-1).lt.0.0d0) 221 > dbl_mb(bmeta(1)+d-1) = rmax 222 end if 223 if (int_mb(indxmeta(1)+s1).eq.2) then 224 if (int_mb(pmeta(1)+d-1).lt.0) 225 > int_mb(pmeta(1)+d-1) = 0 226 if (dbl_mb(ameta(1)+d-1).lt.(-99.0d0)) 227 > dbl_mb(ameta(1)+d-1) = 0.0d0 228 if (dbl_mb(bmeta(1)+d-1).lt.(-99.0d0)) 229 > dbl_mb(bmeta(1)+d-1) = 4.0d0*datan(1.0d0) 230 end if 231 end do 232 end if 233 234* **** load old spline fits **** 235 call dcopy(nxmeta_all,0.0d0,0,dbl_mb(ymeta(1)),1) 236 rtdb_name = 'tamd_ymeta' 237 if (.not.btdb_get(rtdb,rtdb_name,mt_dbl,nxmeta_all, 238 > dbl_mb(ymeta(1)))) then 239 call dcopy(nxmeta_all,0.0d0,0,dbl_mb(ymeta(1)),1) 240 if (boundary(1).gt.0.0d0) then 241 if (taskid.eq.MASTER) 242 > write(luout,'(A,2F11.6)') 243 > " ... setting boundary, w,sigma=", 244 > boundary(1),boundary(2) 245 call tamd_setboundary(boundary(1),boundary(2)) 246 end if 247 if (dabs(sn2(1)).gt.0.0d0) then 248 if (taskid.eq.MASTER) 249 > write(luout,'(A)') " ... setting sn2 surface" 250 call meta_setsn2surface(sn2) 251 end if 252 else 253 if (taskid.eq.MASTER) 254 > write(luout,*) " ... reading tamd data from rtdb" 255 256 !*** re-scale tempered *** 257 fac = 1.0d0 258 if (dTtempered.gt.0.0d0) then 259 if (taskid.eq.MASTER) 260 > write(luout,*) 261 > " ... re-scaling to be consistent with tempering" 262 if (control_Nose()) then 263 temp = control_Nose_Tr() 264 else 265 temp = ion_Temperature() 266 end if 267 fac = dTtempered/(temp+dTtempered) 268 end if 269 do i=1,nxmeta_all 270 dbl_mb(ymeta(1)+i-1) = fac*dbl_mb(ymeta(1)+i-1) 271 end do 272 end if 273 274* **** intialize ztamd and vztamd **** 275 !seed = 99 --- need to have option to read from input 276 !x = util_random(seed) 277 dt = control_time_step() 278 279 !*** read restart ztamd and vztamd from rtdb *** 280 rtdb_name = 'tamd_ztamd_nmeta' 281 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,d)) d = -1 282 ok = .false. 283 if (d.eq.nmeta) then 284 rtdb_name = 'tamd_ztamd' 285 ok = btdb_get(rtdb,rtdb_name,mt_dbl,nmeta,ztamd) 286 rtdb_name = 'tamd_vztamd' 287 ok = ok.and.btdb_get(rtdb,rtdb_name,mt_dbl,nmeta,vztamd) 288 end if 289 !*** read initial ztamd and vztamd from rtdb *** 290 if (.not.ok) then 291 do d=1,nmeta 292 s2 = int_mb(sparammeta(1)+d-1) 293 ztamd(d) = dbl_mb(parammeta(1)+s2+4) 294 vztamd(d) = dbl_mb(parammeta(1)+s2+5) 295 end do 296 end if 297 298 call nwpw_interp_init(nmeta,8) 299 300 if (.not.value) 301 > call errquit('cannot allocate heap memory for tamd',0, 302 > MA_ERR) 303 304 305* **** iniitialize meta_gauss_write **** 306 rtdb_name = 'tamd_gaussian_filename' 307 gauss_filename = 308 > ' ' 309 if(.not.btdb_cget(rtdb,rtdb_name,1,gauss_filename)) 310 > call util_file_prefix('tamd_data',gauss_filename) 311 call tamd_gauss_write_init(gauss_filename) 312 313* **** start nwpw_expression **** 314 call nwpw_expression_start(rtdb) 315 316 317* **** write out header info **** 318 call Parallel_taskid(taskid) 319 value = btdb_parallel(.false.) 320 if (taskid.eq.MASTER) then 321 write(luout,*) 322 write(luout,*) "TAMD parameters:" 323 if (zfixed) 324 > write(luout,'(A)') " - z fixed - determining average force" 325 write(luout,'(A,5x,I6," * inner iterations")') 326 > " - update = ",maxmetacount 327 write(luout,'(A,5x,I6," * inner iterations")') 328 > " - print = ",maxmetaprintcount 329 write(luout,'(A,5x,I6)') " - print_shift = ",metarayshift 330 write(luout,'(A,5x,F6.3)') " - time step = ",dt 331 if (dTtempered.gt.0.0d0) then 332 write(luout,'(A,5x,F11.1)') 333 > " - Tempered dT (K) = ",dTtempered 334 end if 335 do d=1,nmeta 336 s1 = int_mb(sindxmeta(1) +d-1) 337 s2 = int_mb(sparammeta(1)+d-1) 338 ntype = int_mb(indxmeta(1)+s1) 339 gamma = dbl_mb(parammeta(1)+s2) 340 k = dbl_mb(parammeta(1)+s2+1) 341 mass = dbl_mb(parammeta(1)+s2+2) 342 ztemp = dbl_mb(parammeta(1)+s2+3) 343 zinitial(1) = dbl_mb(parammeta(1)+s2+4) 344 zinitial(2) = dbl_mb(parammeta(1)+s2+5) 345 346 if (ntype.eq.1) then 347 write(luout,'(A,5x,2I4, 348 > /6x,A,F13.6,4x,A,F13.6, 349 > /6x,A,F13.6,4x,A,F13.6, 350 > /6x,A,F13.6,F13.6, 351 > /6x,A,F13.6,F13.6,4x,A,I5,4x,A,I2)') 352 > " - Bond Parameters = ", 353 > (int_mb(indxmeta(1)+s1+j),j=1,2), 354 > 'gamma= ',gamma,'k= ',k, 355 > 'mass= ',mass, 356 > 'ztemp=',ztemp, 357 > 'zinitial=',zinitial, 358 > 'range= ',dbl_mb(ameta(1)+d-1),dbl_mb(bmeta(1)+d-1), 359 > 'nrange= ',int_mb(nxmeta(1)+d-1), 360 > 'periodic=',int_mb(pmeta(1)+d-1) 361 362 else if (ntype.eq.2) then 363 write(luout,'(A,3I4, 364 > /6x,A,F13.6,4x,A,F13.6, 365 > /6x,A,F13.6,4x,A,F13.6, 366 > /6x,A,F13.6,F13.6, 367 > /6x,A,F13.6,F13.6,4x,A,I5,4x,A,I2)') 368 > " - Angle Parameters = ", 369 > (int_mb(indxmeta(1)+s1+j),j=1,3), 370 > 'gamma=',gamma,'k=',k, 371 > 'mass= ',mass, 372 > 'ztemp=',ztemp, 373 > 'zinitial=',zinitial, 374 > 'range= ',dbl_mb(ameta(1)+d-1),dbl_mb(bmeta(1)+d-1), 375 > 'nrange= ',int_mb(nxmeta(1)+d-1), 376 > 'periodic=',int_mb(pmeta(1)+d-1) 377 378 else if (ntype.eq.3) then 379 write(luout,'(A,4I4, 380 > /6x,A,F13.6,4x,A,F13.6, 381 > /6x,A,F13.6,4x,A,F13.6, 382 > /6x,A,F13.6,F13.6, 383 > /6x,A,F13.6,F13.6,4x,A,I5,4x,A,I2)') 384 > " - Angle Parameters = ", 385 > (int_mb(indxmeta(1)+s1+j),j=1,4), 386 > 'gamma=',gamma,'k=',k, 387 > 'mass= ',mass, 388 > 'ztemp=',ztemp, 389 > 'zinitial=',zinitial, 390 > 'range= ',dbl_mb(ameta(1)+d-1),dbl_mb(bmeta(1)+d-1), 391 > 'nrange= ',int_mb(nxmeta(1)+d-1), 392 > 'periodic=',int_mb(pmeta(1)+d-1) 393 394 else if (ntype.eq.4) then 395 n1 = int_mb(indxmeta(1)+s1+1) 396 n2 = int_mb(indxmeta(1)+s1+2) 397 n = dbl_mb(parammeta(1)+s2+8) 398 m = dbl_mb(parammeta(1)+s2+9) 399 r0 = dbl_mb(parammeta(1)+s2+10) 400 n12= dbl_mb(parammeta(1)+s2+12) 401 write(luout,1001) (int_mb(indxmeta(1)+s1+3+j-1),j=1,n1) 402 write(luout,1002) (int_mb(indxmeta(1)+s1+3+n1+j-1),j=1,n2) 403 404 write(luout,'(A,/6x,A,F13.6,4x,A,F13.6,4x,A,F13.6, 405 > /6x,A,F13.6, 406 > /6x,A,F13.6,4x,A,F13.6, 407 > /6x,A,F13.6,4x,A,F13.6, 408 > /6x,A,F13.6,F13.6, 409 > /6x,A,F13.6,F13.6, 410 > /6x,A,I5, 411 > /6x,A,I2)') 412 > " - Coordination Number Parameters: ", 413 > 'n= ',n,'m= ',m,'n12= ',n12, 414 > 'r0=',r0, 415 > 'gamma=',gamma,'k=',k, 416 > 'mass= ',mass, 417 > 'ztemp=',ztemp, 418 > 'zinitial=',zinitial, 419 > 'range= ',dbl_mb(ameta(1)+d-1),dbl_mb(bmeta(1)+d-1), 420 > 'ngrid= ',int_mb(nxmeta(1)+d-1), 421 > 'periodic= ',int_mb(pmeta(1)+d-1) 422 if (dbl_mb(parammeta(1)+s2+11).lt.0) then 423 write(luout,'(6x,A)') '- LJ function form' 424 else 425 write(luout,'(6x,A)') '- Sprik function form' 426 end if 427 428 else if (ntype.eq.5) then 429 n1 = int_mb(indxmeta(1)+s1+2) 430 write(luout,1001) int_mb(indxmeta(1)+s1+1) 431 write(luout,1002) (int_mb(indxmeta(1)+s1+3+j-1),j=1,n1) 432 write(luout,'(A,5x, 433 > /6x,A,F13.6,F13.6,F13.6, 434 > /6x,A,F13.6,4x,A,F13.6, 435 > /6x,A,F13.6,4x,A,F13.6, 436 > /6x,A,F13.6,F13.6, 437 > /6x,A,F13.6,F13.6,4x,A,I5,4x,A,I2)') 438 > " - n-plane Parameters = ", 439 > 'normal=', 440 > dbl_mb(parammeta(1)+s2+6), 441 > dbl_mb(parammeta(1)+s2+7), 442 > dbl_mb(parammeta(1)+s2+8), 443 > 'gamma=',gamma,'k=',k, 444 > 'mass= ',dbl_mb(parammeta(1)+s2+2), 445 > 'ztemp=',dbl_mb(parammeta(1)+s2+3), 446 > 'zinitial=',dbl_mb(parammeta(1)+s2+4), 447 > dbl_mb(parammeta(1)+s2+5), 448 > 'range= ',dbl_mb(ameta(1)+d-1),dbl_mb(bmeta(1)+d-1), 449 > 'nrange= ',int_mb(nxmeta(1)+d-1), 450 > 'periodic= ',int_mb(pmeta(1)+d-1) 451 else if (ntype.eq.6) then 452 453 if (int_mb(indxmeta(1)+s1+1).eq.1) then 454 write(luout,'(A,5x,I4, 455 > /6x,A,F13.6,4x,A,F13.6, 456 > /6x,A,F13.6,4x,A,F13.6, 457 > /6x,A,F13.6,F13.6, 458 > /6x,A,F13.6,F13.6,4x,A,I5,4x,A,I2)') 459 > " - X Parameters, Ion= ", 460 > int_mb(indxmeta(1)+s1+2), 461 > 'gamma= ',gamma,'k=',k, 462 > 'mass= ',dbl_mb(parammeta(1)+s2+2), 463 > 'ztemp=',dbl_mb(parammeta(1)+s2+3), 464 > 'zinitial=',dbl_mb(parammeta(1)+s2+4), 465 > dbl_mb(parammeta(1)+s2+5), 466 > 'range= ',dbl_mb(ameta(1)+d-1),dbl_mb(bmeta(1)+d-1), 467 > 'nrange= ',int_mb(nxmeta(1)+d-1), 468 > 'periodic= ',int_mb(pmeta(1)+d-1) 469 else if (int_mb(indxmeta(1)+s1+1).eq.2) then 470 write(luout,'(A,5x,I4, 471 > /6x,A,F13.6,4x,A,F13.6, 472 > /6x,A,F13.6,4x,A,F13.6, 473 > /6x,A,F13.6,F13.6, 474 > /6x,A,F13.6,F13.6,4x,A,I5,4x,A,I2)') 475 > " - Y Parameters, Ion = ", 476 > int_mb(indxmeta(1)+s1+2), 477 > 'gamma= ',gamma,'k=',k, 478 > 'mass= ',dbl_mb(parammeta(1)+s2+2), 479 > 'ztemp=',dbl_mb(parammeta(1)+s2+3), 480 > 'zinitial=',dbl_mb(parammeta(1)+s2+4), 481 > dbl_mb(parammeta(1)+s2+5), 482 > 'range= ',dbl_mb(ameta(1)+d-1),dbl_mb(bmeta(1)+d-1), 483 > 'nrange= ',int_mb(nxmeta(1)+d-1), 484 > 'periodic= ',int_mb(pmeta(1)+d-1) 485 else 486 write(luout,'(A,5x,I4, 487 > /6x,A,F13.6,4x,A,F13.6, 488 > /6x,A,F13.6,4x,A,F13.6, 489 > /6x,A,F13.6,F13.6, 490 > /6x,A,F13.6,F13.6,4x,A,I5,4x,A,I2)') 491 > " - Z Parameters, Ion = ",int_mb(indxmeta(1)+s1+2), 492 > 'gamma= ',gamma,'k=',k, 493 > 'mass= ',dbl_mb(parammeta(1)+s2+2), 494 > 'ztemp=',dbl_mb(parammeta(1)+s2+3), 495 > 'zinitial=',dbl_mb(parammeta(1)+s2+4), 496 > dbl_mb(parammeta(1)+s2+5), 497 > 'range= ',dbl_mb(ameta(1)+d-1),dbl_mb(bmeta(1)+d-1), 498 > 'nrange= ',int_mb(nxmeta(1)+d-1), 499 > 'periodic= ',int_mb(pmeta(1)+d-1) 500 end if 501 502 else if (ntype.eq.7) then 503 write(luout,'(A,5x,A,I4,4x,A,I4, 504 > /6x,A,F13.6,4x,A,F13.6, 505 > /6x,A,F13.6,4x,A,F13.6, 506 > /6x,A,F13.6,F13.6, 507 > /6x,A,F13.6,F13.6,4x,A,I5,4x,A,I2)') 508 > " - local_density_trace Parameters = ", 509 > 'atom index=',int_mb(indxmeta(1)+s1+1), 510 > 'l=',int_mb(indxmeta(1)+s1+2), 511 > 'gamma=',gamma,'k=',k, 512 > 'mass= ',dbl_mb(parammeta(1)+s2+2), 513 > 'ztemp=',dbl_mb(parammeta(1)+s2+3), 514 > 'zinitial=',dbl_mb(parammeta(1)+s2+4), 515 > dbl_mb(parammeta(1)+s2+5) 516c > 'range= ',dbl_mb(ameta(1)+d-1),dbl_mb(bmeta(1)+d-1), 517c > 'nrange= ',int_mb(nxmeta(1)+d-1), 518c > 'periodic=',int_mb(pmeta(1)+d-1) 519 520 else if (ntype.eq.8) then 521 write(luout,'(A, 522 > /6x,A,2I4, 523 > /6x,A,2I4, 524 > /6x,A,2F13.6, 525 > /6x,A,2F13.6, 526 > /6x,A,F13.6,4x,A,F13.6, 527 > /6x,A,F13.6,4x,A,F13.6, 528 > /6x,A,F13.6,F13.6, 529 > /6x,A,F13.6,F13.6,4x,A,I5,4x,A,I2)') 530 > " - 2Bonds Parameters = ", 531 > 'bond 1 indexes=',(int_mb(indxmeta(1)+s1+j),j=1,2), 532 > 'bond 2 indexes=',(int_mb(indxmeta(1)+s1+j),j=3,4), 533 > 'bond center 1 =',(dbl_mb(parammeta(1)+s2+j),j=2,3), 534 > 'bond center 2 =',(dbl_mb(parammeta(1)+s2+j),j=4,5), 535 > 'gamma=',gamma,'k=',k, 536 > 'mass= ',dbl_mb(parammeta(1)+s2+2), 537 > 'ztemp=',dbl_mb(parammeta(1)+s2+3), 538 > 'zinitial=',dbl_mb(parammeta(1)+s2+4), 539 > dbl_mb(parammeta(1)+s2+5) 540c > 'range= ',dbl_mb(ameta(1)+d-1),dbl_mb(bmeta(1)+d-1), 541c > 'nrange= ',int_mb(nxmeta(1)+d-1), 542c > 'periodic=',int_mb(pmeta(1)+d-1) 543 544 else if (ntype.eq.9) then 545 write(luout,'(A,5x,A,I4,4x,A,I4,4x,A,I4, 546 > /6x,A,F13.6,4x,A,F13.6, 547 > /6x,A,F13.6,4x,A,F13.6, 548 > /6x,A,F13.6,F13.6, 549 > /6x,A,F13.6,F13.6,4x,A,I5,4x,A,I2)') 550 > " - local_density_trace2_diff Parameters = ", 551 > 'atom index1=',int_mb(indxmeta(1)+s1+1), 552 > 'atom index2=',int_mb(indxmeta(1)+s1+2), 553 > 'l=',int_mb(indxmeta(1)+s1+3), 554 > 'gamma=',gamma,'k=',k, 555 > 'mass= ',dbl_mb(parammeta(1)+s2+2), 556 > 'ztemp=',dbl_mb(parammeta(1)+s2+3), 557 > 'zinitial=',dbl_mb(parammeta(1)+s2+4), 558 > dbl_mb(parammeta(1)+s2+5) 559c > 'range= ',dbl_mb(ameta(1)+d-1),dbl_mb(bmeta(1)+d-1), 560c > 'nrange= ',int_mb(nxmeta(1)+d-1), 561c > 'periodic=',int_mb(pmeta(1)+d-1) 562 563 else if (ntype.eq.10) then 564 write(luout,'(A, 565 > /6x,A,3I4, 566 > /6x,A,F13.6,4x,A,F13.6, 567 > /6x,A,F13.6,4x,A,F13.6, 568 > /6x,A,F13.6,F13.6, 569 > /6x,A,F13.6,F13.6,4x,A,I5,4x,A,I2)') 570 > " - Bond_Difference2 Parameters = ", 571 > 'bond indexes=',(int_mb(indxmeta(1)+s1+j),j=1,3), 572 > 'gamma=',gamma,'k=',k, 573 > 'mass= ',dbl_mb(parammeta(1)+s2+2), 574 > 'ztemp=',dbl_mb(parammeta(1)+s2+3), 575 > 'zinitial=',dbl_mb(parammeta(1)+s2+4), 576 > dbl_mb(parammeta(1)+s2+5) 577c > 'range= ',dbl_mb(ameta(1)+d-1),dbl_mb(bmeta(1)+d-1), 578c > 'nrange= ',int_mb(nxmeta(1)+d-1), 579c > 'periodic=',int_mb(pmeta(1)+d-1) 580 581 else if (ntype.eq.12) then 582 call nwpw_expression_eqnstring(rtdb, 583 > int_mb(indxmeta(1)+s1+1), 584 > eqnstring) 585 write(luout,'(A,/6x,A,A, 586 > /6x,A,F13.6,4x,A,F13.6, 587 > /6x,A,F13.6,4x,A,F13.6, 588 > /6x,A,F13.6,F13.6, 589 > /6x,A,F13.6,F13.6,4x,A,I5,4x,A,I2)') 590 > ' - Equation Parameters: ', 591 > 'equation= ',trim(eqnstring), 592 > 'gamma=',gamma,'k=',k, 593 > 'mass= ',dbl_mb(parammeta(1)+s2+2), 594 > 'ztemp=',dbl_mb(parammeta(1)+s2+3), 595 > 'zinitial=',dbl_mb(parammeta(1)+s2+4), 596 > dbl_mb(parammeta(1)+s2+5) 597c > 'range= ',dbl_mb(ameta(1)+d-1),dbl_mb(bmeta(1)+d-1), 598c > 'nrange= ',int_mb(nxmeta(1)+d-1), 599c > 'periodic=',int_mb(pmeta(1)+d-1) 600 601 602 else if (ntype.eq.13) then 603 n1 = int_mb(indxmeta(1)+s1+1) 604 n2 = int_mb(indxmeta(1)+s1+2) 605 n = dbl_mb(parammeta(1)+s2+6) 606 write(luout,1005) (int_mb(indxmeta(1)+s1+3+j-1),j=1,n1) 607 write(luout,1006) (int_mb(indxmeta(1)+s1+3+n1+j-1),j=1,n2) 608 write(luout,'(A,/6x,A,F13.6, 609 > /6x,A,F13.6,4x,A,F13.6, 610 > /6x,A,F13.6,4x,A,F13.6, 611 > /6x,A,F13.6,F13.6, 612 > /6x,A,F13.6,F13.6,4x,A,I5,4x,A,I2)') 613 > " - Group_Distance Parameters: ", 614 > 'n= ',n, 615 > 'gamma=',gamma,'k=',k, 616 > 'mass= ',dbl_mb(parammeta(1)+s2+2), 617 > 'ztemp=',dbl_mb(parammeta(1)+s2+3), 618 > 'zinitial=',dbl_mb(parammeta(1)+s2+4), 619 > dbl_mb(parammeta(1)+s2+5), 620 > 'range= ',dbl_mb(ameta(1)+d-1),dbl_mb(bmeta(1)+d-1), 621 > 'nrange= ',int_mb(nxmeta(1)+d-1), 622 > 'periodic=',int_mb(pmeta(1)+d-1) 623 624 else if (ntype.eq.14) then 625 n1 = int_mb(indxmeta(1)+s1+1) 626 write(luout,'(A, 627 > /6x,A,F13.6,4x,A,F13.6, 628 > /6x,A,F13.6,4x,A,F13.6, 629 > /6x,A,F13.6,F13.6, 630 > /6x,A,F13.6,F13.6,4x,A,I5,4x,A,I2)') 631 > " - Bondings Parameters: ", 632 > 'gamma=',gamma,'k=',k, 633 > 'mass= ',dbl_mb(parammeta(1)+s2+2), 634 > 'ztemp=',dbl_mb(parammeta(1)+s2+3), 635 > 'zinitial=',dbl_mb(parammeta(1)+s2+4), 636 > dbl_mb(parammeta(1)+s2+5), 637 > 'range= ',dbl_mb(ameta(1)+d-1),dbl_mb(bmeta(1)+d-1), 638 > 'nrange= ',int_mb(nxmeta(1)+d-1), 639 > 'periodic=',int_mb(pmeta(1)+d-1) 640 641 write(luout,'(A)')" coefficient index1 index2 :" 642 do j=1,n1 643 write(luout,'(F12.6,2I7)') dbl_mb(parammeta(1)+s2+5+j), 644 > int_mb(indxmeta(1)+s1+2*j), 645 > int_mb(indxmeta(1)+s1+2*j+1) 646 end do 647 648 end if 649 650 end do 651 652 1001 FORMAT(3x,"- Coorination Number (Index1) :",10I5) 653 1002 FORMAT(3x,"- Coorination Number (Index2) :",10I5) 654 1005 FORMAT(2x,"- Group_Distance Indexes (Index1) :",10I5) 655 1006 FORMAT(2x,"- Group_Distance Indexes (Index2) :",10I5) 656 657 end if 658 value = btdb_parallel(.true.) 659 660 end if 661 return 662 end 663 664c ********************************************** 665c * * 666c * tamd_finalize * 667c * * 668c ********************************************** 669 subroutine tamd_finalize(rtdb) 670 implicit none 671 integer rtdb 672 673#include "bafdecls.fh" 674#include "btdb.fh" 675#include "errquit.fh" 676#include "tamd.fh" 677 678* **** local variables **** 679 integer taskid,MASTER 680 parameter (MASTER=0) 681 682 logical value 683 integer i,j,yshift,k 684 real*8 temp,fac 685 character*80 filename,rtdb_name 686 character*60 gauss_filename 687 character*255 full_filename 688 689* **** external functions **** 690 logical control_Nose 691 external control_Nose 692 character*7 c_index_name 693 external c_index_name 694 real*8 control_Nose_Tr,ion_Temperature 695 external control_Nose_Tr,ion_Temperature 696 697 call Parallel_taskid(taskid) 698 699 if (tamdfound) then 700 701* **** print out potentials ****** 702 !metaraycount = metaraycount + metarayshift 703 !call tamd_print_potentials(0) 704 705 value=.true. 706 rtdb_name='tamd_metacount' 707 value= value.and.btdb_put(rtdb,rtdb_name,mt_int,1,metacount) 708 rtdb_name='tamd_metaraycount' 709 value= value.and.btdb_put(rtdb,rtdb_name,mt_int,1,metaraycount) 710 711 if (nmeta.gt.0) then 712 713 !*** un-scale tempered *** 714 fac = 1.0d0 715 if (dTtempered.gt.0.0d0) then 716 if (control_Nose()) then 717 temp = control_Nose_Tr() 718 else 719 temp = ion_Temperature() 720 end if 721 fac = (temp+dTtempered)/dTtempered 722 end if 723 do i=1,nxmeta_all 724 dbl_mb(ymeta(1)+i-1) = fac*dbl_mb(ymeta(1)+i-1) 725 end do 726 727* **** write out current ztamd and vztamd **** 728 rtdb_name = 'tamd_ztamd_nmeta' 729 value = value.and. 730 > btdb_put(rtdb,rtdb_name,mt_int,1,nmeta) 731 rtdb_name = 'tamd_ztamd' 732 value = value.and. 733 > btdb_put(rtdb,rtdb_name,mt_dbl,nmeta,ztamd) 734 rtdb_name = 'tamd_vztamd' 735 value = value.and. 736 > btdb_put(rtdb,rtdb_name,mt_dbl,nmeta,vztamd) 737 738* **** write out current spline fits **** 739 rtdb_name = 'tamd_ymeta' 740 value = value.and.btdb_put(rtdb,rtdb_name,mt_dbl, 741 > nxmeta_all,dbl_mb(ymeta(1))) 742 743 value = value.and.BA_free_heap(pmeta(2)) 744 value = value.and.BA_free_heap(ameta(2)) 745 value = value.and.BA_free_heap(bmeta(2)) 746 value = value.and.BA_free_heap(nxmeta(2)) 747 value = value.and.BA_free_heap(indxmeta(2)) 748 value = value.and.BA_free_heap(parammeta(2)) 749 value = value.and.BA_free_heap(sindxmeta(2)) 750 value = value.and.BA_free_heap(sparammeta(2)) 751 value = value.and.BA_free_heap(ymeta(2)) 752 if (.not.value) 753 > call errquit('tamd_end: cannot free heap',0,MA_ERR) 754 755 call nwpw_interp_end() 756 call tamd_gauss_write_end() 757 call nwpw_expression_end(rtdb) 758 759 760 761* **** analyze tamd_data_file **** 762 if (zfixed) then 763 rtdb_name = 'tamd_gaussian_filename' 764 gauss_filename = 765 > ' ' 766 if(.not.btdb_cget(rtdb,rtdb_name,1,gauss_filename)) 767 > call util_file_prefix('tamd_data',gauss_filename) 768 call tamd_gauss_analyze(rtdb,gauss_filename) 769 end if 770 771 end if 772 773 end if 774 775 return 776 end 777 778 779c ********************************************** 780c * * 781c * tamd_collective * 782c * * 783c ********************************************** 784 real*8 function tamd_collective(indx,param,ispin,ne,psi1) 785 implicit none 786 integer indx(*) 787 real*8 param(*) 788 integer ispin,ne(2) 789 real*8 psi1(*) 790 791#include "bafdecls.fh" 792 793* **** local variables **** 794 integer i,j,ntype,n1,n2,nion 795 real*8 x1,y1,z1,x2,y2,z2,x3,y3,z3,dx,dy,dz,r,r1,r2,r3 796 real*8 dx1,dy1,dz1,dx3,dy3,dz3,theta,n,m,r0,f,la,lb 797 real*8 nx,ny,nz,x4,y4,z4,dx2,dy2,dz2,d1a,d2a,d1b,d2b,n12 798 799* **** external functions **** 800 integer ion_rion_ptr,ion_nion 801 external ion_rion_ptr,ion_nion 802 real*8 ion_rion,psp_ld_trace 803 external ion_rion,psp_ld_trace 804 real*8 nwpw_expression_f 805 external nwpw_expression_f 806 807 808 ntype = indx(1) 809 810 f = 0.0d0 811c **** bond distance *** 812 if (ntype.eq.1) then 813 x1 = ion_rion(1,indx(2)) 814 y1 = ion_rion(2,indx(2)) 815 z1 = ion_rion(3,indx(2)) 816 x2 = ion_rion(1,indx(3)) 817 y2 = ion_rion(2,indx(3)) 818 z2 = ion_rion(3,indx(3)) 819 dx = x1-x2 820 dy = y1-y2 821 dz = z1-z2 822 call lattice_min_difference(dx,dy,dz) 823 r = dsqrt(dx**2 + dy**2 + dz**2) 824 f = r 825 826c **** bond angle *** 827 else if (ntype.eq.2) then 828 x1 = ion_rion(1,indx(2)) 829 y1 = ion_rion(2,indx(2)) 830 z1 = ion_rion(3,indx(2)) 831 x2 = ion_rion(1,indx(3)) 832 y2 = ion_rion(2,indx(3)) 833 z2 = ion_rion(3,indx(3)) 834 x3 = ion_rion(1,indx(4)) 835 y3 = ion_rion(2,indx(4)) 836 z3 = ion_rion(3,indx(4)) 837 dx1 = x1-x2 838 dy1 = y1-y2 839 dz1 = z1-z2 840 call lattice_min_difference(dx1,dy1,dz1) 841 r1 = dsqrt(dx1**2 + dy1**2 + dz1**2) 842 dx3 = x3-x2 843 dy3 = y3-y2 844 dz3 = z3-z2 845 call lattice_min_difference(dx3,dy3,dz3) 846 r3 = dsqrt(dx3**2 + dy3**2 + dz3**2) 847 theta = (dx1*dx3 + dy1*dy3 + dz1*dz3)/(r1*r3) 848 if (theta.gt.1.0d0) theta = 1.0d0 849 if (theta.lt.-1.0d0) theta = -1.0d0 850 theta = dacos(theta) 851 f = theta 852 853c **** bond dihedral **** 854 else if (ntype.eq.3) then 855 x1 = ion_rion(1,indx(2)) 856 y1 = ion_rion(2,indx(2)) 857 z1 = ion_rion(3,indx(2)) 858 x2 = ion_rion(1,indx(3)) 859 y2 = ion_rion(2,indx(3)) 860 z2 = ion_rion(3,indx(3)) 861 x3 = ion_rion(1,indx(4)) 862 y3 = ion_rion(2,indx(4)) 863 z3 = ion_rion(3,indx(4)) 864 x4 = ion_rion(1,indx(5)) 865 y4 = ion_rion(2,indx(5)) 866 z4 = ion_rion(3,indx(5)) 867 dx1 = x2-x1 868 dy1 = y2-y1 869 dz1 = z2-z1 870 call lattice_min_difference(dx1,dy1,dz1) 871 dx2 = x3-x2 872 dy2 = y3-y2 873 dz2 = z3-z2 874 call lattice_min_difference(dx2,dy2,dz2) 875 dx3 = x4-x3 876 dy3 = y4-y3 877 dz3 = z4-z3 878 call lattice_min_difference(dx3,dy3,dz3) 879 880 dy = dx1*(dy2*dz3-dy3*dz2) 881 > + dy1*(dz2*dx3-dz3*dx2) 882 > + dz1*(dx2*dy3-dx3*dy2) 883 dy = dy*dsqrt(dx2**2 + dy2**2 + dz2**2) 884 885 dx = (dy2*dz3 - dy3*dz2)*(dy1*dz2 - dy2*dz1) 886 > + (dz2*dx3 - dz3*dx2)*(dz1*dx2 - dz2*dx1) 887 > + (dx2*dy3 - dx3*dy2)*(dx1*dy2 - dx2*dy1) 888 f = datan2(dy,dx) 889 890c **** coordination number **** 891 else if (ntype.eq.4) then 892 n1=indx(2) 893 n2=indx(3) 894 if (param(12).lt.0.0d0) then 895 f = 0.0d0 896 n = param(9) 897 m = param(10) 898 r0 = param(11) 899 n12= param(13) 900 do i=1,n1 901 do j=1,n2 902 x1 = ion_rion(1,indx(3+i)) 903 y1 = ion_rion(2,indx(3+i)) 904 z1 = ion_rion(3,indx(3+i)) 905 x2 = ion_rion(1,indx(3+n1+j)) 906 y2 = ion_rion(2,indx(3+n1+j)) 907 z2 = ion_rion(3,indx(3+n1+j)) 908 dx1 = x1-x2 909 dy1 = y1-y2 910 dz1 = z1-z2 911 call lattice_min_difference(dx1,dy1,dz1) 912 r1 = dsqrt(dx1**2 + dy1**2 + dz1**2) 913 f = f + (1.0d0-(r1/r0)**n)/(1.0d0-(r1/r0)**m) 914 end do 915 end do 916 f = f/n12 917 else 918 f = 0.0d0 919 n = param(9) 920 m = param(10) 921 r0 = param(11) 922 n12= param(13) 923 do i=1,n1 924 do j=1,n2 925 x1 = ion_rion(1,indx(3+i)) 926 y1 = ion_rion(2,indx(3+i)) 927 z1 = ion_rion(3,indx(3+i)) 928 x2 = ion_rion(1,indx(3+n1+j)) 929 y2 = ion_rion(2,indx(3+n1+j)) 930 z2 = ion_rion(3,indx(3+n1+j)) 931 dx1 = x1-x2 932 dy1 = y1-y2 933 dz1 = z1-z2 934 call lattice_min_difference(dx1,dy1,dz1) 935 r1 = dsqrt(dx1**2 + dy1**2 + dz1**2) 936c f = f + (1.0d0-(r1/r0)**n)/(1.0d0-(r1/r0)**m) 937 f = f + 1.d0/(1.d0+dexp(n*(r1-r0))) 938 end do 939 end do 940 f = f/n12 941 end if 942 943c **** n-plane **** 944 else if (ntype.eq.5) then 945 f = 0.0d0 946 nx = param(7) 947 ny = param(8) 948 nz = param(9) 949 x1 = ion_rion(1,indx(2)) 950 y1 = ion_rion(2,indx(2)) 951 z1 = ion_rion(3,indx(2)) 952 n1=indx(3) 953 do i=1,n1 954 x2 = ion_rion(1,indx(3+i)) 955 y2 = ion_rion(2,indx(3+i)) 956 z2 = ion_rion(3,indx(3+i)) 957 dx = x1-x2 958 dy = y1-y2 959 dz = z1-z2 960 call lattice_min_difference(dx,dy,dz) 961 r = nx*dx + ny*dy + nz*dz 962 f = f + dsqrt(r*r) 963 end do 964 f = f/dble(n1) 965 966* **** x,y,z **** 967 else if (ntype.eq.6) then 968 if (indx(2).eq.1) then 969 f = ion_rion(1,indx(3)) 970 else if (indx(2).eq.2) then 971 f = ion_rion(2,indx(3)) 972 else 973 f = ion_rion(3,indx(3)) 974 end if 975 976* **** ld_trace **** 977 else if (ntype.eq.7) then 978 f = psp_ld_trace(indx(2),indx(3),ispin,ne,psi1) 979 980* *** 2bonds **** 981 else if (ntype.eq.8) then 982 x1 = ion_rion(1,indx(2)) 983 y1 = ion_rion(2,indx(2)) 984 z1 = ion_rion(3,indx(2)) 985 x2 = ion_rion(1,indx(3)) 986 y2 = ion_rion(2,indx(3)) 987 z2 = ion_rion(3,indx(3)) 988 x3 = ion_rion(1,indx(4)) 989 y3 = ion_rion(2,indx(4)) 990 z3 = ion_rion(3,indx(4)) 991 x4 = ion_rion(1,indx(5)) 992 y4 = ion_rion(2,indx(5)) 993 z4 = ion_rion(3,indx(5)) 994 dx1 = x2-x1 995 dy1 = y2-y1 996 dz1 = z2-z1 997 call lattice_min_difference(dx1,dy1,dz1) 998 dx2 = x4-x3 999 dy2 = y4-y3 1000 dz2 = z4-z3 1001 call lattice_min_difference(dx2,dy2,dz2) 1002 d1a = param(9) 1003 d2a = param(10) 1004 d1b = param(11) 1005 d2b = param(12) 1006 r1 = dsqrt(dx1**2 + dy1**2 + dz1**2) 1007 r2 = dsqrt(dx2**2 + dy2**2 + dz2**2) 1008 la = dsqrt((r1-d1a)**2 + (r2-d2a)**2) 1009 lb = dsqrt((r1-d1b)**2 + (r2-d2b)**2) 1010 f = la/(la+lb) 1011 1012* **** ld_trace2diff **** 1013 else if (ntype.eq.9) then 1014 f = psp_ld_trace(indx(2),indx(4),ispin,ne,psi1) 1015 > - psp_ld_trace(indx(3),indx(4),ispin,ne,psi1) 1016 1017* *** bond_difference2 **** 1018 else if (ntype.eq.10) then 1019 x1 = ion_rion(1,indx(2)) 1020 y1 = ion_rion(2,indx(2)) 1021 z1 = ion_rion(3,indx(2)) 1022 x2 = ion_rion(1,indx(3)) 1023 y2 = ion_rion(2,indx(3)) 1024 z2 = ion_rion(3,indx(3)) 1025 x3 = ion_rion(1,indx(4)) 1026 y3 = ion_rion(2,indx(4)) 1027 z3 = ion_rion(3,indx(4)) 1028 dx1 = x2-x1 1029 dy1 = y2-y1 1030 dz1 = z2-z1 1031 call lattice_min_difference(dx1,dy1,dz1) 1032 r1 = (dx1**2 + dy1**2 + dz1**2) 1033 dx2 = x3-x1 1034 dy2 = y3-y1 1035 dz2 = z3-z1 1036 call lattice_min_difference(dx2,dy2,dz2) 1037 r2 = (dx2**2 + dy2**2 + dz2**2) 1038 f = r1 - r2 1039 1040 else if (ntype.eq.12) then 1041 nion = ion_nion() 1042 f = nwpw_expression_f(indx(2),nion,dbl_mb(ion_rion_ptr())) 1043 1044c **** group distance **** 1045 else if (ntype.eq.13) then 1046 n1=indx(2) 1047 n2=indx(3) 1048 f = 0.0d0 1049 n = param(7) 1050 n12 = dble(n1*n2) 1051 do i=1,n1 1052 do j=1,n2 1053 x1 = ion_rion(1,indx(3+i)) 1054 y1 = ion_rion(2,indx(3+i)) 1055 z1 = ion_rion(3,indx(3+i)) 1056 x2 = ion_rion(1,indx(3+n1+j)) 1057 y2 = ion_rion(2,indx(3+n1+j)) 1058 z2 = ion_rion(3,indx(3+n1+j)) 1059 dx1 = x1-x2 1060 dy1 = y1-y2 1061 dz1 = z1-z2 1062 call lattice_min_difference(dx1,dy1,dz1) 1063 r1 = dsqrt(dx1**2 + dy1**2 + dz1**2) 1064 f = f + (1.0/r1)**n 1065 end do 1066 end do 1067 f = (f/n12)**(-1.0d0/n) 1068 1069c **** bondings **** 1070 else if (ntype.eq.14) then 1071 n1=indx(2) 1072 f = 0.0d0 1073 do i=1,n1 1074 x1 = ion_rion(1,indx(2*i+1)) 1075 y1 = ion_rion(2,indx(2*i+1)) 1076 z1 = ion_rion(3,indx(2*i+1)) 1077 x2 = ion_rion(1,indx(2*i+2)) 1078 y2 = ion_rion(2,indx(2*i+2)) 1079 z2 = ion_rion(3,indx(2*i+2)) 1080 dx1 = x1-x2 1081 dy1 = y1-y2 1082 dz1 = z1-z2 1083 call lattice_min_difference(dx1,dy1,dz1) 1084 r1 = dsqrt(dx1**2 + dy1**2 + dz1**2) 1085 f = f + param(i+6)*r1 1086 end do 1087 end if 1088 1089 tamd_collective = f 1090 return 1091 end 1092 1093c ********************************************** 1094c * * 1095c * tamd_collective_force * 1096c * * 1097c ********************************************** 1098 subroutine tamd_collective_force(indx,param,dv, 1099 > ispin,ne,psi1,fmeta_psi1, 1100 > move,fmeta) 1101 implicit none 1102 integer indx(*) 1103 real*8 param(*),dv 1104 integer ispin,ne(2) 1105 real*8 psi1(*),fmeta_psi1(*) 1106 logical move 1107 real*8 fmeta(3,*) 1108 1109#include "bafdecls.fh" 1110 1111* **** local variables **** 1112 integer i,j,ntype,n1,n2,nion 1113 real*8 x1,y1,z1,r1,vx1,vx2,vy1,vy2,vz1,vz2 1114 real*8 x2,y2,z2,aa,a11,a12,a22,r,df,n,m,r0,rn,rm 1115 real*8 x3,y3,z3,r3,ctheta,stheta,denom 1116 real*8 dx,dy,dz,dx1,dy1,dz1,dx3,dy3,dz3,nx,ny,nz,la,lb 1117 real*8 x4,y4,z4,dx2,dy2,dz2,r2,df1,df2,d1a,d2a,d1b,d2b 1118 real*8 dfla,dflb,dlar1,dlar2,dlbr1,dlbr2,n12,f 1119 1120* **** external functions **** 1121 integer ion_rion_ptr,ion_nion 1122 external ion_rion_ptr,ion_nion 1123 real*8 ion_rion 1124 external ion_rion 1125 1126 ntype = indx(1) 1127 1128 if (move) then 1129 1130c **** bond distance *** 1131 if (ntype.eq.1) then 1132 x1 = ion_rion(1,indx(2)) 1133 y1 = ion_rion(2,indx(2)) 1134 z1 = ion_rion(3,indx(2)) 1135 x2 = ion_rion(1,indx(3)) 1136 y2 = ion_rion(2,indx(3)) 1137 z2 = ion_rion(3,indx(3)) 1138 dx = x1-x2 1139 dy = y1-y2 1140 dz = z1-z2 1141 call lattice_min_difference(dx,dy,dz) 1142 r = dsqrt(dx**2 + dy**2 + dz**2) 1143 1144 fmeta(1,indx(2)) = fmeta(1,indx(2)) - (dx/r)*dv 1145 fmeta(2,indx(2)) = fmeta(2,indx(2)) - (dy/r)*dv 1146 fmeta(3,indx(2)) = fmeta(3,indx(2)) - (dz/r)*dv 1147 fmeta(1,indx(3)) = fmeta(1,indx(3)) + (dx/r)*dv 1148 fmeta(2,indx(3)) = fmeta(2,indx(3)) + (dy/r)*dv 1149 fmeta(3,indx(3)) = fmeta(3,indx(3)) + (dz/r)*dv 1150 1151c **** bond angle *** 1152 else if (ntype.eq.2) then 1153 x1 = ion_rion(1,indx(2)) 1154 y1 = ion_rion(2,indx(2)) 1155 z1 = ion_rion(3,indx(2)) 1156 x2 = ion_rion(1,indx(3)) 1157 y2 = ion_rion(2,indx(3)) 1158 z2 = ion_rion(3,indx(3)) 1159 x3 = ion_rion(1,indx(4)) 1160 y3 = ion_rion(2,indx(4)) 1161 z3 = ion_rion(3,indx(4)) 1162 dx1 = x1-x2 1163 dy1 = y1-y2 1164 dz1 = z1-z2 1165 call lattice_min_difference(dx1,dy1,dz1) 1166 r1 = dsqrt(dx1**2 + dy1**2 + dz1**2) 1167 dx3 = x3-x2 1168 dy3 = y3-y2 1169 dz3 = z3-z2 1170 call lattice_min_difference(dx3,dy3,dz3) 1171 r3 = dsqrt(dx3**2 + dy3**2 + dz3**2) 1172 denom = r1*r3 1173 if (denom.gt.1.0d-11) then 1174 ctheta = (dx1*dx3 + dy1*dy3 + dz1*dz3)/(denom) 1175 if (ctheta.gt.1.0d0) ctheta = 1.0d0 1176 if (ctheta.lt.-1.0d0) ctheta = -1.0d0 1177 stheta = dsqrt(1.0d0-ctheta*ctheta) 1178 if (stheta.lt.0.001d0) stheta = 0.001d0 1179 stheta = 1.0d0/stheta 1180 1181 aa = dv*stheta 1182 a11 = aa*ctheta/r1 1183 a12 = -aa/(denom) 1184 a22 = aa*ctheta/r3 1185 1186 vx1 = a11*dx1 + a12*dx3 1187 vx2 = a22*dx3 + a12*dx1 1188 1189 vy1 = a11*dy1 + a12*dy3 1190 vy2 = a22*dy3 + a12*dy1 1191 1192 vz1 = a11*dz1 + a12*dz3 1193 vz2 = a22*dz3 + a12*dz1 1194 1195 fmeta(1,indx(2)) = fmeta(1,indx(2)) - vx1 1196 fmeta(2,indx(2)) = fmeta(2,indx(2)) - vy1 1197 fmeta(3,indx(2)) = fmeta(3,indx(2)) - vz1 1198 1199 fmeta(1,indx(3)) = fmeta(1,indx(3)) + vx1 + vx2 1200 fmeta(2,indx(3)) = fmeta(2,indx(3)) + vy1 + vy2 1201 fmeta(3,indx(3)) = fmeta(3,indx(3)) + vz1 + vz2 1202 1203 fmeta(1,indx(4)) = fmeta(1,indx(4)) - vx2 1204 fmeta(2,indx(4)) = fmeta(2,indx(4)) - vy2 1205 fmeta(3,indx(4)) = fmeta(3,indx(4)) - vz2 1206 end if 1207 1208c **** bond dihedral **** 1209 else if (ntype.eq.3) then 1210 x1 = ion_rion(1,indx(2)) 1211 y1 = ion_rion(2,indx(2)) 1212 z1 = ion_rion(3,indx(2)) 1213 x2 = ion_rion(1,indx(3)) 1214 y2 = ion_rion(2,indx(3)) 1215 z2 = ion_rion(3,indx(3)) 1216 x3 = ion_rion(1,indx(4)) 1217 y3 = ion_rion(2,indx(4)) 1218 z3 = ion_rion(3,indx(4)) 1219 x4 = ion_rion(1,indx(5)) 1220 y4 = ion_rion(2,indx(5)) 1221 z4 = ion_rion(3,indx(5)) 1222 dx1 = x2-x1 1223 dy1 = y2-y1 1224 dz1 = z2-z1 1225 call lattice_min_difference(dx1,dy1,dz1) 1226 dx2 = x3-x2 1227 dy2 = y3-y2 1228 dz2 = z3-z2 1229 call lattice_min_difference(dx2,dy2,dz2) 1230 dx3 = x4-x3 1231 dy3 = y4-y3 1232 dz3 = z4-z3 1233 call lattice_min_difference(dx3,dy3,dz3) 1234 1235 dy = dx1*(dy2*dz3-dy3*dz2) 1236 > + dy1*(dz2*dx3-dz3*dx2) 1237 > + dz1*(dx2*dy3-dx3*dy2) 1238 dy = dy*dsqrt(dx2**2 + dy2**2 + dz2**2) 1239 1240 dx = (dy2*dz3 - dy3*dz2)*(dy1*dz2 - dy2*dz1) 1241 > + (dz2*dx3 - dz3*dx2)*(dz1*dx2 - dz2*dx1) 1242 > + (dx2*dy3 - dx3*dy2)*(dx1*dy2 - dx2*dy1) 1243 r2 = datan2(dy,dx) 1244 vx1 = -dy/(dx**2 + dy**2) 1245 vy1 = dx/(dx**2 + dy**2) 1246 1247c dphi/dx (dx/db1 db1/dr1 + dx/db2 db2/dr1 + dx/db3 db3/dr1) 1248c dphi/dy (dy/db1 db1/dr1 + dy/db2 db2/dr1 + dy/db3 db3/dr1) 1249 1250 fmeta(1,indx(2)) = fmeta(1,indx(2)) 1251 fmeta(2,indx(2)) = fmeta(2,indx(2)) 1252 fmeta(3,indx(2)) = fmeta(3,indx(2)) 1253 1254 fmeta(1,indx(3)) = fmeta(1,indx(3)) 1255 fmeta(2,indx(3)) = fmeta(2,indx(3)) 1256 fmeta(3,indx(3)) = fmeta(3,indx(3)) 1257 1258 fmeta(1,indx(4)) = fmeta(1,indx(4)) 1259 fmeta(2,indx(4)) = fmeta(2,indx(4)) 1260 fmeta(3,indx(4)) = fmeta(3,indx(4)) 1261 1262 fmeta(1,indx(5)) = fmeta(1,indx(5)) 1263 fmeta(2,indx(5)) = fmeta(2,indx(5)) 1264 fmeta(3,indx(5)) = fmeta(3,indx(5)) 1265 1266 1267 1268c **** coordination number *** 1269 else if (ntype.eq.4) then 1270 n1 = indx(2) 1271 n2 = indx(3) 1272 if (param(12).lt.0.0d0) then 1273 n = param(9) 1274 m = param(10) 1275 r0 = param(11) 1276 n12= param(13) 1277 do i=1,n1 1278 do j=1,n2 1279 x1 = ion_rion(1,indx(3+i)) 1280 y1 = ion_rion(2,indx(3+i)) 1281 z1 = ion_rion(3,indx(3+i)) 1282 x2 = ion_rion(1,indx(3+n1+j)) 1283 y2 = ion_rion(2,indx(3+n1+j)) 1284 z2 = ion_rion(3,indx(3+n1+j)) 1285 dx = x1-x2 1286 dy = y1-y2 1287 dz = z1-z2 1288 call lattice_min_difference(dx,dy,dz) 1289 r = dsqrt(dx**2 + dy**2 + dz**2) 1290 rn = (1.0d0-(r/r0)**n) 1291 rm = (1.0d0-(r/r0)**m) 1292 df = (-n*rm/r0*(r/r0)**(n-1) + m*rn/r0*(r/r0)**(m-1)) 1293 > / (rm**2) 1294 df = df/n12 1295 fmeta(1,indx(3+i)) = fmeta(1,indx(3+i)) - (dx/r)*df*dv 1296 fmeta(2,indx(3+i)) = fmeta(2,indx(3+i)) - (dy/r)*df*dv 1297 fmeta(3,indx(3+i)) = fmeta(3,indx(3+i)) - (dz/r)*df*dv 1298 fmeta(1,indx(3+n1+j)) = fmeta(1,indx(3+n1+j)) 1299 > + (dx/r)*df*dv 1300 fmeta(2,indx(3+n1+j)) = fmeta(2,indx(3+n1+j)) 1301 > + (dy/r)*df*dv 1302 fmeta(3,indx(3+n1+j)) = fmeta(3,indx(3+n1+j)) 1303 > + (dz/r)*df*dv 1304 end do 1305 end do 1306 else 1307 n = param(9) 1308 m = param(10) 1309 r0 = param(11) 1310 n12= param(13) 1311 do i=1,n1 1312 do j=1,n2 1313 x1 = ion_rion(1,indx(3+i)) 1314 y1 = ion_rion(2,indx(3+i)) 1315 z1 = ion_rion(3,indx(3+i)) 1316 x2 = ion_rion(1,indx(3+n1+j)) 1317 y2 = ion_rion(2,indx(3+n1+j)) 1318 z2 = ion_rion(3,indx(3+n1+j)) 1319 dx = x1-x2 1320 dy = y1-y2 1321 dz = z1-z2 1322 call lattice_min_difference(dx,dy,dz) 1323 r = dsqrt(dx**2 + dy**2 + dz**2) 1324c rn = (1.0d0-(r/r0)**n) 1325c rm = (1.0d0-(r/r0)**m) 1326c df = (-n*rm/r0*(r/r0)**(n-1) + m*rn/r0*(r/r0)**(m-1)) 1327c > / (rm**2) 1328 rn = 1.d0+dexp(n*(r-r0)) 1329 df = -n*(rn-1.d0)/(rn*rn) 1330 df = df/n12 1331 fmeta(1,indx(3+i)) = fmeta(1,indx(3+i)) - (dx/r)*df*dv 1332 fmeta(2,indx(3+i)) = fmeta(2,indx(3+i)) - (dy/r)*df*dv 1333 fmeta(3,indx(3+i)) = fmeta(3,indx(3+i)) - (dz/r)*df*dv 1334 fmeta(1,indx(3+n1+j)) = fmeta(1,indx(3+n1+j)) 1335 > + (dx/r)*df*dv 1336 fmeta(2,indx(3+n1+j)) = fmeta(2,indx(3+n1+j)) 1337 > + (dy/r)*df*dv 1338 fmeta(3,indx(3+n1+j)) = fmeta(3,indx(3+n1+j)) 1339 > + (dz/r)*df*dv 1340 end do 1341 end do 1342 endif 1343 1344c **** n-plane **** 1345 else if (ntype.eq.5) then 1346 nx = param(9) 1347 ny = param(10) 1348 nz = param(11) 1349 x1 = ion_rion(1,indx(2)) 1350 y1 = ion_rion(2,indx(2)) 1351 z1 = ion_rion(3,indx(2)) 1352 n1=indx(3) 1353 do i=1,n1 1354 x2 = ion_rion(1,indx(3+i)) 1355 y2 = ion_rion(2,indx(3+i)) 1356 z2 = ion_rion(3,indx(3+i)) 1357 dx = x1-x2 1358 dy = y1-y2 1359 dz = z1-z2 1360 call lattice_min_difference(dx,dy,dz) 1361 r = nx*dx + ny*dy + nz*dz 1362 r0 = dsqrt(r*r) 1363 df = (r/r0)*dv/dble(n1) 1364 fmeta(1,indx(2)) = fmeta(1,indx(2)) - nx*df 1365 fmeta(2,indx(2)) = fmeta(2,indx(2)) - ny*df 1366 fmeta(3,indx(2)) = fmeta(3,indx(2)) - nz*df 1367 1368 fmeta(1,indx(3+i)) = fmeta(1,indx(3+i)) + nx*df 1369 fmeta(2,indx(3+i)) = fmeta(2,indx(3+i)) + ny*df 1370 fmeta(3,indx(3+i)) = fmeta(3,indx(3+i)) + nz*df 1371 end do 1372 1373* **** x,y,z **** 1374 else if (ntype.eq.6) then 1375 if (indx(2).eq.1) then 1376 fmeta(1,indx(3)) = fmeta(1,indx(3)) - dv 1377 else if (indx(2).eq.2) then 1378 fmeta(2,indx(3)) = fmeta(2,indx(3)) - dv 1379 else 1380 fmeta(3,indx(3)) = fmeta(3,indx(3)) - dv 1381 end if 1382 1383c **** ld_trace *** 1384 else if (ntype.eq.7) then 1385 call psp_ld_trace_gradient(indx(2),indx(3),ispin,ne,psi1, 1386 > .false.,param, 1387 > (-dv),fmeta_psi1, 1388 > move,fmeta) 1389 1390c **** 2bonds **** 1391 else if (ntype.eq.8) then 1392 x1 = ion_rion(1,indx(2)) 1393 y1 = ion_rion(2,indx(2)) 1394 z1 = ion_rion(3,indx(2)) 1395 x2 = ion_rion(1,indx(3)) 1396 y2 = ion_rion(2,indx(3)) 1397 z2 = ion_rion(3,indx(3)) 1398 x3 = ion_rion(1,indx(4)) 1399 y3 = ion_rion(2,indx(4)) 1400 z3 = ion_rion(3,indx(4)) 1401 x4 = ion_rion(1,indx(5)) 1402 y4 = ion_rion(2,indx(5)) 1403 z4 = ion_rion(3,indx(5)) 1404 dx1 = x2-x1 1405 dy1 = y2-y1 1406 dz1 = z2-z1 1407 call lattice_min_difference(dx1,dy1,dz1) 1408 r1 = dsqrt(dx1**2 + dy1**2 + dz1**2) 1409 dx2 = x4-x3 1410 dy2 = y4-y3 1411 dz2 = z4-z3 1412 call lattice_min_difference(dx2,dy2,dz2) 1413 r2 = dsqrt(dx2**2 + dy2**2 + dz2**2) 1414 d1a = param(9) 1415 d2a = param(10) 1416 d1b = param(11) 1417 d2b = param(12) 1418 la = dsqrt((r1-d1a)**2 + (r2-d2a)**2) 1419 lb = dsqrt((r1-d1b)**2 + (r2-d2b)**2) 1420 if (dabs(la).lt.1.0d-6) then 1421 dlar1 = dtanh((r1-d1a)/1.0d-7) 1422 dlar2 = dtanh((r2-d2a)/1.0d-7) 1423 else 1424 dlar1 = (r1-d1a)/la 1425 dlar2 = (r2-d2a)/la 1426 end if 1427 if (dabs(lb).lt.1.0d-6) then 1428 dlbr1 = dtanh((r1-d1b)/1.0d-7) 1429 dlbr2 = dtanh((r2-d2b)/1.0d-7) 1430 else 1431 dlbr1 = (r1-d1b)/lb 1432 dlbr2 = (r2-d2b)/lb 1433 end if 1434 !f = la/(la+lb) 1435 dfla = 1.0d0/(la+lb) - la/(la+lb)**2 1436 dflb = -la/(la+lb)**2 1437 df1 = dfla*dlar1 + dflb*dlbr1 1438 df2 = dfla*dlar2 + dflb*dlbr2 1439 1440 fmeta(1,indx(2)) = fmeta(1,indx(2)) + (dx1/r1)*dv*df1 1441 fmeta(2,indx(2)) = fmeta(2,indx(2)) + (dy1/r1)*dv*df1 1442 fmeta(3,indx(2)) = fmeta(3,indx(2)) + (dz1/r1)*dv*df1 1443 1444 fmeta(1,indx(3)) = fmeta(1,indx(3)) - (dx1/r1)*dv*df1 1445 fmeta(2,indx(3)) = fmeta(2,indx(3)) - (dy1/r1)*dv*df1 1446 fmeta(3,indx(3)) = fmeta(3,indx(3)) - (dz1/r1)*dv*df1 1447 1448 fmeta(1,indx(4)) = fmeta(1,indx(4)) + (dx2/r2)*dv*df2 1449 fmeta(2,indx(4)) = fmeta(2,indx(4)) + (dy2/r2)*dv*df2 1450 fmeta(3,indx(4)) = fmeta(3,indx(4)) + (dz2/r2)*dv*df2 1451 1452 fmeta(1,indx(5)) = fmeta(1,indx(5)) - (dx2/r2)*dv*df2 1453 fmeta(2,indx(5)) = fmeta(2,indx(5)) - (dy2/r2)*dv*df2 1454 fmeta(3,indx(5)) = fmeta(3,indx(5)) - (dz2/r2)*dv*df2 1455 1456c **** ld_trace2diff *** 1457 else if (ntype.eq.9) then 1458 call psp_ld_trace_gradient(indx(2),indx(4),ispin,ne,psi1, 1459 > .false.,param, 1460 > (-dv),fmeta_psi1, 1461 > move,fmeta) 1462 call psp_ld_trace_gradient(indx(3),indx(4),ispin,ne,psi1, 1463 > .false.,param, 1464 > (dv),fmeta_psi1, 1465 > move,fmeta) 1466 1467c **** bond_difference2 **** 1468 else if (ntype.eq.10) then 1469 x1 = ion_rion(1,indx(2)) 1470 y1 = ion_rion(2,indx(2)) 1471 z1 = ion_rion(3,indx(2)) 1472 x2 = ion_rion(1,indx(3)) 1473 y2 = ion_rion(2,indx(3)) 1474 z2 = ion_rion(3,indx(3)) 1475 x3 = ion_rion(1,indx(4)) 1476 y3 = ion_rion(2,indx(4)) 1477 z3 = ion_rion(3,indx(4)) 1478 dx1 = x2-x1 1479 dy1 = y2-y1 1480 dz1 = z2-z1 1481 call lattice_min_difference(dx1,dy1,dz1) 1482 r1 = (dx1**2 + dy1**2 + dz1**2) 1483 dx2 = x3-x1 1484 dy2 = y3-y1 1485 dz2 = z3-z1 1486 call lattice_min_difference(dx2,dy2,dz2) 1487 r2 = (dx2**2 + dy2**2 + dz2**2) 1488 !f = r1 - r2 1489 df1 = 1 1490 df2 = -1 1491 fmeta(1,indx(2)) = fmeta(1,indx(2)) +2.0d0*(dx1*df1+dx2*df2)*dv 1492 fmeta(2,indx(2)) = fmeta(2,indx(2)) +2.0d0*(dy1*df1+dy2*df2)*dv 1493 fmeta(3,indx(2)) = fmeta(3,indx(2)) +2.0d0*(dz1*df1+dz2*df2)*dv 1494 1495 fmeta(1,indx(3)) = fmeta(1,indx(3)) - (2.0d0*dx1*df1)*dv 1496 fmeta(2,indx(3)) = fmeta(2,indx(3)) - (2.0d0*dy1*df1)*dv 1497 fmeta(3,indx(3)) = fmeta(3,indx(3)) - (2.0d0*dz1*df1)*dv 1498 1499 fmeta(1,indx(4)) = fmeta(1,indx(4)) - (2.0d0*dx2*df2)*dv 1500 fmeta(2,indx(4)) = fmeta(2,indx(4)) - (2.0d0*dy2*df2)*dv 1501 fmeta(3,indx(4)) = fmeta(3,indx(4)) - (2.0d0*dz2*df2)*dv 1502 1503 else if (ntype.eq.12) then 1504 nion = ion_nion() 1505 call nwpw_expression_fion(dv,indx(2), 1506 > nion,dbl_mb(ion_rion_ptr()), 1507 > fmeta) 1508 1509c **** group distances **** 1510 else if (ntype.eq.13) then 1511 n1=indx(2) 1512 n2=indx(3) 1513 f = 0.0d0 1514 n = param(7) 1515 n12 = dble(n1*n2) 1516 do i=1,n1 1517 do j=1,n2 1518 x1 = ion_rion(1,indx(3+i)) 1519 y1 = ion_rion(2,indx(3+i)) 1520 z1 = ion_rion(3,indx(3+i)) 1521 x2 = ion_rion(1,indx(3+n1+j)) 1522 y2 = ion_rion(2,indx(3+n1+j)) 1523 z2 = ion_rion(3,indx(3+n1+j)) 1524 dx1 = x1-x2 1525 dy1 = y1-y2 1526 dz1 = z1-z2 1527 call lattice_min_difference(dx1,dy1,dz1) 1528 r1 = dsqrt(dx1**2 + dy1**2 + dz1**2) 1529 f = f + (1.0/r1)**n 1530 end do 1531 end do 1532 f = (f/n12)**(-1.0d0/n) 1533 1534 do i=1,n1 1535 do j=1,n2 1536 x1 = ion_rion(1,indx(3+i)) 1537 y1 = ion_rion(2,indx(3+i)) 1538 z1 = ion_rion(3,indx(3+i)) 1539 x2 = ion_rion(1,indx(3+n1+j)) 1540 y2 = ion_rion(2,indx(3+n1+j)) 1541 z2 = ion_rion(3,indx(3+n1+j)) 1542 dx = x1-x2 1543 dy = y1-y2 1544 dz = z1-z2 1545 call lattice_min_difference(dx,dy,dz) 1546 r1 = dsqrt(dx**2 + dy**2 + dz**2) 1547 df = ((f/r1)**(n+1))/n12 1548 1549 fmeta(1,indx(3+i)) = fmeta(1,indx(3+i)) - (dx/r1)*df*dv 1550 fmeta(2,indx(3+i)) = fmeta(2,indx(3+i)) - (dy/r1)*df*dv 1551 fmeta(3,indx(3+i)) = fmeta(3,indx(3+i)) - (dz/r1)*df*dv 1552 fmeta(1,indx(3+n1+j)) = fmeta(1,indx(3+n1+j)) 1553 > + (dx/r1)*df*dv 1554 fmeta(2,indx(3+n1+j)) = fmeta(2,indx(3+n1+j)) 1555 > + (dy/r1)*df*dv 1556 fmeta(3,indx(3+n1+j)) = fmeta(3,indx(3+n1+j)) 1557 > + (dz/r1)*df*dv 1558 1559 1560 end do 1561 end do 1562 1563c **** bondings **** 1564 else if (ntype.eq.14) then 1565 n1=indx(2) 1566 do i=1,n1 1567 x1 = ion_rion(1,indx(2*i+1)) 1568 y1 = ion_rion(2,indx(2*i+1)) 1569 z1 = ion_rion(3,indx(2*i+1)) 1570 x2 = ion_rion(1,indx(2*i+2)) 1571 y2 = ion_rion(2,indx(2*i+2)) 1572 z2 = ion_rion(3,indx(2*i+2)) 1573 dx1 = x1-x2 1574 dy1 = y1-y2 1575 dz1 = z1-z2 1576 call lattice_min_difference(dx1,dy1,dz1) 1577 r1 = dsqrt(dx1**2 + dy1**2 + dz1**2) 1578 fmeta(1,indx(2*i+1)) = fmeta(1,indx(2*i+1)) 1579 > - param(i+6)*dx1/r1*dv 1580 fmeta(2,indx(2*i+1)) = fmeta(2,indx(2*i+1)) 1581 > - param(i+6)*dy1/r1*dv 1582 fmeta(3,indx(2*i+1)) = fmeta(3,indx(2*i+1)) 1583 > - param(i+6)*dz1/r1*dv 1584 fmeta(1,indx(2*i+2)) = fmeta(1,indx(2*i+2)) 1585 > + param(i+6)*dx1/r1*dv 1586 fmeta(2,indx(2*i+2)) = fmeta(2,indx(2*i+2)) 1587 > + param(i+6)*dy1/r1*dv 1588 fmeta(3,indx(2*i+2)) = fmeta(3,indx(2*i+2)) 1589 > + param(i+6)*dz1/r1*dv 1590 end do 1591 1592 1593 end if 1594 1595* ***** move = .false. ****** 1596 else 1597 if (ntype.eq.7) then 1598 call psp_ld_trace_gradient(indx(2),indx(3),ispin,ne,psi1, 1599 > .false.,param, 1600 > (-dv),fmeta_psi1, 1601 > move,fmeta) 1602 else if (ntype.eq.9) then 1603 call psp_ld_trace_gradient(indx(2),indx(4),ispin,ne,psi1, 1604 > .false.,param, 1605 > (-dv),fmeta_psi1, 1606 > move,fmeta) 1607 call psp_ld_trace_gradient(indx(3),indx(4),ispin,ne,psi1, 1608 > .false.,param, 1609 > (dv),fmeta_psi1, 1610 > move,fmeta) 1611 end if 1612 end if 1613 1614 return 1615 end 1616 1617 1618 1619* *********************************************** 1620* * * 1621* * tamd_update * 1622* * * 1623* *********************************************** 1624 subroutine tamd_update(ispin,ne,psi1,E) 1625 implicit none 1626 integer ispin,ne(2) 1627 real*8 psi1(*) 1628 real*8 E(*) 1629 1630#include "bafdecls.fh" 1631#include "errquit.fh" 1632#include "tamd.fh" 1633 1634* **** local variables **** 1635 double precision kb 1636 parameter (kb=3.16679d-6) 1637 1638 integer d,i,j(10),s1,s2 1639 real*8 theta(10),r,x,w,ww,sigma,v,expv 1640 1641* **** external functions **** 1642 real*8 tamd_collective,tamd_energy 1643 external tamd_collective,tamd_energy 1644 1645 if (tamdfound) then 1646 1647 if (nmeta.gt.0) then 1648 do d=1,nmeta 1649 s1 = int_mb(sindxmeta(1) +d-1) 1650 s2 = int_mb(sparammeta(1)+d-1) 1651 theta(d) = tamd_collective(int_mb(indxmeta(1)+s1), 1652 > dbl_mb(parammeta(1)+s2), 1653 > ispin,ne,psi1) 1654 end do 1655 call tamd_gauss_write(nmeta,theta,ztamd,vztamd,E) 1656 end if 1657 1658 end if 1659 1660 return 1661 end 1662 1663 1664* *********************************************** 1665* * * 1666* * tamd_force * 1667* * * 1668* *********************************************** 1669 subroutine tamd_force(ispin,ne,psi1,fmeta_psi1,move,fmeta) 1670 implicit none 1671 integer ispin,ne(2) 1672 real*8 psi1(*),fmeta_psi1(*) 1673 logical move 1674 real*8 fmeta(3,*) 1675 1676#include "bafdecls.fh" 1677#include "errquit.fh" 1678#include "tamd.fh" 1679 1680* **** local variables **** 1681 integer d,s1,s2,p 1682 real*8 r(10),dv(10),k,gam,mz,ztemp,a,b 1683 1684* **** external functions **** 1685 real*8 tamd_collective 1686 external tamd_collective 1687 1688 if (tamdfound) then 1689 1690 if (nmeta.gt.0) then 1691 do d=1,nmeta 1692 s1 = int_mb(sindxmeta(1) +d-1) 1693 s2 = int_mb(sparammeta(1)+d-1) 1694 gam = dbl_mb(parammeta(1)+s2) 1695 k = dbl_mb(parammeta(1)+s2+1) 1696 mz = dbl_mb(parammeta(1)+s2+2) 1697 ztemp = dbl_mb(parammeta(1)+s2+3) 1698 r(d) = tamd_collective(int_mb(indxmeta(1)+s1), 1699 > dbl_mb(parammeta(1)+s2), 1700 > ispin,ne,psi1) 1701 1702 dv(d) = k*(r(d)-ztamd(d)) 1703 call tamd_collective_force(int_mb(indxmeta(1) +s1), 1704 > dbl_mb(parammeta(1)+s2), 1705 > dv(d), 1706 > ispin,ne,psi1,fmeta_psi1, 1707 > move,fmeta) 1708 1709c **** update ztamd and vztamd **** 1710 if (.not.zfixed) then 1711 call tamd_zupdate(dt,gam,k,mz,ztemp,r(d), 1712 > ztamd(d),vztamd(d)) 1713 p = int_mb(pmeta(1)+d-1) 1714 a = dbl_mb(ameta(1)+d-1) 1715 b = dbl_mb(bmeta(1)+d-1) 1716 if (ztamd(d).lt.a) ztamd(d) = a + p*(b-a) 1717 if (ztamd(d).gt.b) ztamd(d) = b + p*(a-b) 1718 end if 1719 end do 1720 1721 end if 1722 1723 1724 end if 1725 1726 return 1727 end 1728 1729* *********************************************** 1730* * * 1731* * tamd_zupdate * 1732* * * 1733* *********************************************** 1734 subroutine tamd_zupdate(dt,gam,k,mz,ztemp,theta,z,vz) 1735 implicit none 1736 real*8 dt,gam,k,mz,ztemp,theta,z,vz 1737 1738 1739* **** local variables **** 1740 real*8 kb 1741 parameter (kb=3.16679d-6) 1742 1743 integer i 1744 real*8 sigma,Az,zf, zf_new,z_new,vz_new,re,rn 1745 real*8 uuu1,uuu2 1746 1747* **** external functions **** 1748 real*8 util_random 1749 external util_random 1750 1751 sigma = dsqrt(2.0d0*kb*ztemp*gam)/mz 1752 zf= 1.0d0*k*(theta-z) 1753 1754 uuu1=util_random(0) 1755 uuu2=util_random(0) 1756 uuu1=2.0d0*3.1415927d0*uuu1 1757 uuu2=dsqrt(-2.0d0*dlog(uuu2)) 1758 re = uuu2*dcos(uuu1) 1759 rn = uuu2*dsin(uuu1) 1760 1761 Az=0.5d0*dt*dt*(zf/mz-gam*vz)+sigma*dt**1.5d0*(0.5d0*re+ 1762 > 0.5d0/dsqrt(3.0d0)*rn) 1763 1764 z_new = z+dt*vz+Az 1765 1766 zf_new = 1.0d0*k*(theta-z_new) 1767 1768 vz_new = vz + 0.5d0*dt*(zf+zf_new)/mz-dt*gam*vz+dt**0.5d0* 1769 > sigma*re-gam*Az 1770 1771 z = z_new 1772 vz = vz_new 1773 return 1774 end 1775 1776 1777* *********************************************** 1778* * * 1779* * tamd_energy * 1780* * * 1781* *********************************************** 1782 real*8 function tamd_energy(ispin,ne,psi1) 1783 implicit none 1784 integer ispin,ne(2) 1785 real*8 psi1(*) 1786 1787#include "bafdecls.fh" 1788#include "errquit.fh" 1789#include "tamd.fh" 1790 1791* **** local variables **** 1792 integer d,s1,s2 1793 real*8 r(10),v,k 1794 1795* **** external functions **** 1796 real*8 nwpw_interp,tamd_collective 1797 external nwpw_interp,tamd_collective 1798 1799 v = 0.0d0 1800 if (tamdfound) then 1801 1802 if (nmeta.gt.0) then 1803 do d=1,nmeta 1804 s1 = int_mb(sindxmeta(1) +d-1) 1805 s2 = int_mb(sparammeta(1)+d-1) 1806 k = dbl_mb(parammeta(1)+s2+1) 1807 r(d) = tamd_collective(int_mb(indxmeta(1)+s1), 1808 > dbl_mb(parammeta(1)+s2), 1809 > ispin,ne,psi1) 1810 v = v + 0.5d0*k*(r(d)-ztamd(d))**2 1811 end do 1812 call Parallel_Brdcst_value(0,v) 1813 1814 end if 1815 1816 end if 1817 1818 tamd_energy = v 1819 return 1820 end 1821 1822 1823 1824 1825 1826* *********************************************** 1827* * * 1828* * tamd_energypotential * 1829* * * 1830* *********************************************** 1831 subroutine tamd_energypotential(ispin,ne,psi1,e,p) 1832 implicit none 1833 integer ispin,ne(2) 1834 real*8 psi1(*) 1835 real*8 e,p 1836 1837#include "bafdecls.fh" 1838#include "errquit.fh" 1839#include "tamd.fh" 1840 1841* **** local variables **** 1842 integer d,s1,s2,ntype 1843 real*8 r(10),dv(10),v,k 1844 1845* **** external functions **** 1846 real*8 nwpw_interp,tamd_collective 1847 external nwpw_interp,tamd_collective 1848 1849 e = 0.0d0 1850 p = 0.0d0 1851 if (tamdfound) then 1852 1853 if (nmeta.gt.0) then 1854 do d=1,nmeta 1855 s1 = int_mb(sindxmeta(1) +d-1) 1856 s2 = int_mb(sparammeta(1)+d-1) 1857 k = dbl_mb(parammeta(1)+s2+1) 1858 r(d) = tamd_collective(int_mb(indxmeta(1)+s1), 1859 > dbl_mb(parammeta(1)+s2), 1860 > ispin,ne,psi1) 1861 e = e + 0.5d0*k*(r(d)-ztamd(d))**2 1862 end do 1863 call Parallel_Brdcst_value(0,e) 1864 1865 end if 1866 1867 end if 1868 1869 return 1870 end 1871 1872 1873* *********************************************** 1874* * * 1875* * tamd_found * 1876* * * 1877* *********************************************** 1878 logical function tamd_found() 1879 implicit none 1880 1881#include "tamd.fh" 1882 1883 tamd_found = tamdfound 1884 return 1885 end 1886 1887 1888c ********************************************** 1889c * * 1890c * tamd_print_potentials * 1891c * * 1892c ********************************************** 1893 subroutine tamd_print_potentials(icount) 1894 implicit none 1895 integer icount 1896 1897#include "bafdecls.fh" 1898#include "errquit.fh" 1899#include "tamd.fh" 1900 1901* **** local variables **** 1902 integer taskid,MASTER 1903 parameter (MASTER=0) 1904 1905 integer d,i,j(10) 1906 real*8 x(10),r,v,dv(10),fac,temp 1907 character*80 filename 1908 character*255 full_filename 1909 1910* **** external functions **** 1911 logical control_Nose 1912 external control_Nose 1913 character*7 c_index_name 1914 external c_index_name 1915 real*8 nwpw_interp,control_Nose_Tr,ion_Temperature 1916 external nwpw_interp,control_Nose_Tr,ion_Temperature 1917 1918 call Parallel_taskid(taskid) 1919 1920 if (tamdfound) then 1921 1922* **** determine tempered tamd factor **** 1923 fac = 1.0d0 1924 if (dTtempered.gt.0.0d0) then 1925 if (control_Nose()) then 1926 temp = control_Nose_Tr() 1927 else 1928 temp = ion_Temperature() 1929 end if 1930 fac = (temp+dTtempered)/dTtempered 1931 end if 1932 1933 1934* **** print out potentials ****** 1935 if (nmeta.gt.0) then 1936 if (taskid.eq.MASTER) then 1937 if (icount.gt.0) then 1938 filename = "meta"//c_index_name(icount)//".dat" 1939 else 1940 filename = "meta"//".dat" 1941 end if 1942 1943 call util_file_name_noprefix(filename,.false., 1944 > .false., 1945 > full_filename) 1946 open(unit=53,file=full_filename,form='formatted') 1947 1948 write(53,'(A,A)') "#tamd filename=",filename 1949 write(53,'(A,I6)') "#nmeta ",nmeta 1950 write(53,'(A,10I6)') "#nxmeta ", 1951 > (int_mb(nxmeta(1)+i-1),i=1,nmeta) 1952 write(53,'(A,10I6)') "#sindxmeta ", 1953 > (int_mb(sindxmeta(1)+i-1),i=1,nmeta) 1954 write(53,'(A,10I6)') "#sparammeta ", 1955 > (int_mb(sparammeta(1)+i-1),i=1,nmeta) 1956 write(53,'(A,10F18.12)') "#ameta ", 1957 > (dbl_mb(ameta(1)+i-1),i=1,nmeta) 1958 write(53,'(A,10F18.12)') "#bmeta ", 1959 > (dbl_mb(bmeta(1)+i-1),i=1,nmeta) 1960 write(53,'(A,10I3)') "#pmeta ", 1961 > (int_mb(pmeta(1)+i-1),i=1,nmeta) 1962 write(53,'(A,I8)') "#nindxmeta ",nindxmeta 1963 write(53,'(A,500I6)') "#indxmeta ", 1964 > (int_mb(indxmeta(1)+i-1),i=1,nindxmeta) 1965 write(53,'(A,I8)') "#nparammeta ",nparammeta 1966 write(53,'(A,500F18.12)') "#parammeta ", 1967 > (dbl_mb(parammeta(1)+i-1),i=1,nparammeta) 1968 write(53,'(A,I8)') "#maxmetacount ",maxmetacount 1969 write(53,'(A,I8)') "#metacount ",metacount 1970 write(53,'(A,I8)') "#metarayshift ",metarayshift 1971 write(53,'(A,I8)') "#metaraycount ",metaraycount 1972 write(53,'(A,I8)')"#maxmetaprintcount ",maxmetaprintcount 1973 write(53,'(A,F18.12)') "#dTempered ",dTtempered 1974 write(53,'(A,I8)') "#nxmeta_all ",nxmeta_all 1975 1976 do d=1,nmeta 1977 j(d) = 0 1978 end do 1979 do i=1,nxmeta_all 1980 do d=1,nmeta 1981 x(d) = dbl_mb(ameta(1)+d-1) 1982 > + j(d)*(dbl_mb(bmeta(1)+d-1)-dbl_mb(ameta(1)+d-1)) 1983 > /dble(int_mb(nxmeta(1)+d-1)-1) 1984 end do 1985 write(53,'(12F15.6)') 1986 > (x(d),d=1,nmeta),dbl_mb(ymeta(1)+i-1)*fac 1987 j(1) = j(1) + 1 1988 if (j(1).ge.int_mb(nxmeta(1))) then 1989 write(53,*) 1990 end if 1991 do d=1,nmeta-1 1992 if (j(d).ge.int_mb(nxmeta(1)+d-1)) then 1993 j(d) = 0 1994 j(d+1) = j(d+1)+1 1995 end if 1996 end do 1997 end do 1998 1999 close(53) 2000 2001 2002c filename = "test.dat" 2003c call util_file_name_noprefix(filename,.false., 2004c > .false., 2005c > full_filename) 2006c open(unit=53,file=full_filename,form='formatted') 2007c do i=1,501 2008c x(1) = 1.0d0 + (i-1)*8.0d0/501.0d0 2009c 2010c v = nwpw_interp(nmeta,int_mb(nxmeta(1)), 2011c > dbl_mb(ameta(1)),dbl_mb(bmeta(1)), 2012c > dbl_mb(ymeta(1)),x) 2013c call nwpw_dinterp(nmeta,int_mb(nxmeta(1)), 2014c > dbl_mb(ameta(1)),dbl_mb(bmeta(1)), 2015c > dbl_mb(ymeta(1)),x,dv) 2016c 2017c write(53,'(12F15.6)') x(1),v,dv(1) 2018c 2019c end do 2020c close(53) 2021 2022 end if 2023 end if 2024 2025 end if 2026 2027 return 2028 end 2029 2030 2031 2032 2033* *********************************************** 2034* * * 2035* * tamd_setboundary * 2036* * * 2037* *********************************************** 2038 subroutine tamd_setboundary(bheight,bsigma) 2039 implicit none 2040 real*8 bheight,bsigma 2041 2042#include "bafdecls.fh" 2043#include "errquit.fh" 2044#include "tamd.fh" 2045 2046* **** local variables **** 2047 integer d,i,j(10) 2048 real*8 r0,r,x,w,ww,wa,wb,ra,rb 2049 2050 if (tamdfound) then 2051 if (nmeta.gt.0) then 2052 do d=1,nmeta 2053 j(d) = 0 2054 end do 2055 do i=1,nxmeta_all 2056 ww = 0.0d0 2057 do d=1,nmeta 2058 if (pmeta(d).eq.0) then 2059 x = dbl_mb(ameta(1)+d-1) 2060 > + j(d)*(dbl_mb(bmeta(1)+d-1)-dbl_mb(ameta(1)+d-1)) 2061 > /dble(int_mb(nxmeta(1)+d-1)-1) 2062 else 2063 x = dbl_mb(ameta(1)+d-1) 2064 > + j(d)*(dbl_mb(bmeta(1)+d-1)-dbl_mb(ameta(1)+d-1)) 2065 > /dble(int_mb(nxmeta(1)+d-1)) 2066 end if 2067 ra= dbl_mb(ameta(1)+d-1) 2068 rb= dbl_mb(bmeta(1)+d-1) 2069 wa = dexp(-0.5d0*((x-ra)/bsigma)**2) 2070 wb = dexp(-0.5d0*((x-rb)/bsigma)**2) 2071 ww = ww + wa + wb 2072 end do 2073 2074 dbl_mb(ymeta(1)+i-1) = dbl_mb(ymeta(1)+i-1) + bheight*ww 2075 2076 j(1) = j(1) + 1 2077 do d=1,nmeta-1 2078 if (j(d).ge.int_mb(nxmeta(1)+d-1)) then 2079 j(d) = 0 2080 j(d+1) = j(d+1)+1 2081 end if 2082 end do 2083 end do 2084 end if 2085 end if 2086 return 2087 end 2088 2089 2090 2091* *********************************************** 2092* * * 2093* * tamd_gauss_write_init * 2094* * * 2095* *********************************************** 2096 subroutine tamd_gauss_write_init(filename) 2097 implicit none 2098 character*(*) filename 2099 2100#include "bafdecls.fh" 2101#include "stdio.fh" 2102#include "tamd.fh" 2103 2104 integer MASTER 2105 parameter (MASTER=0) 2106 2107 logical found,found_bak 2108 integer taskid,l1,l2,i 2109 character*255 full_filename,full_bak 2110 2111 call Parallel_taskid(taskid) 2112 2113 2114* **** produce TAMD_DATA FILE **** 2115 if (taskid.eq.MASTER) then 2116 2117 call util_file_name_noprefix(filename,.false., 2118 > .false., 2119 > full_filename) 2120 2121* **** check for backup file *** 2122 call util_file_name_noprefix('TAMD-Z99-bak',.false., 2123 > .false., 2124 > full_bak) 2125 inquire(file=full_bak,exist=found_bak) 2126 if (found_bak) then 2127 write(luout,*) 2128 write(luout,*) "TAMD-Z99-bak exists:" 2129 l1=index(full_bak,' ') 2130 l2=index(full_filename,' ') 2131 write(luout,*) " Copying ",full_bak(1:l2), 2132 > " to ",full_filename(1:l2) 2133 write(*,*) 2134 call util_file_copy(full_bak,full_filename) 2135 end if 2136 inquire(file=full_filename,exist=found) 2137 if (found) then 2138 2139* **** make a new backup file *** 2140 call util_file_copy(full_filename,full_bak) 2141 2142 open(unit=54,file=full_filename,form='formatted', 2143 > status='old') 2144 do while (found) 2145 read(54,*,end=100) 2146 end do 2147 100 continue 2148#if defined(FUJITSU_SOLARIS) || defined(SOLARIS) || defined(__crayx1) || defined(GCC46) 2149 backspace 54 2150#endif 2151 else 2152 open(unit=54,file=full_filename,form='formatted', 2153 > status='new') 2154 2155 write(54,'(A,A40)') 2156 > "#tamd data filename=",filename 2157 write(54,'(A,I6)') "#nmeta ",nmeta 2158 write(54,'(A,10I6)') "#nxmeta ", 2159 > (int_mb(nxmeta(1)+i-1),i=1,nmeta) 2160 write(54,'(A,10I6)') "#sindxmeta ", 2161 > (int_mb(sindxmeta(1)+i-1),i=1,nmeta) 2162 write(54,'(A,10I6)') "#sparammeta ", 2163 > (int_mb(sparammeta(1)+i-1),i=1,nmeta) 2164 write(54,'(A,10F18.12)') "#ameta ", 2165 > (dbl_mb(ameta(1)+i-1),i=1,nmeta) 2166 write(54,'(A,10F18.12)') "#bmeta ", 2167 > (dbl_mb(bmeta(1)+i-1),i=1,nmeta) 2168 write(54,'(A,10I3)') "#pmeta ", 2169 > (int_mb(pmeta(1)+i-1),i=1,nmeta) 2170 write(54,'(A,I8)') "#nindxmeta ",nindxmeta 2171 write(54,'(A,500I6)') "#indxmeta ", 2172 > (int_mb(indxmeta(1)+i-1),i=1,nindxmeta) 2173 write(54,'(A,I8)') "#nparammeta ",nparammeta 2174 write(54,'(A,500F18.12)') "#parammeta ", 2175 > (dbl_mb(parammeta(1)+i-1),i=1,nparammeta) 2176 write(54,'(A,I8)') "#maxmetacount ",maxmetacount 2177 write(54,'(A,I8)') "#metacount ",metacount 2178 write(54,'(A,I8)') "#metarayshift ",metarayshift 2179 write(54,'(A,I8)') "#metaraycount ",metaraycount 2180 write(54,'(A,I8)') "#maxmetaprintcount ",maxmetaprintcount 2181 write(54,'(A,A)') 2182 > "# (theta(i),i=1,n),(z(i),i=1,n),(vz(i),i=1,n),", 2183 > "E(2)-E(33)+E(34),E(3),E(4),E(33),E(34)" 2184 end if 2185 end if 2186 2187 return 2188 end 2189 2190 2191 2192* *********************************************** 2193* * * 2194* * tamd_gauss_write_end * 2195* * * 2196* *********************************************** 2197 subroutine tamd_gauss_write_end() 2198 implicit none 2199 2200 integer MASTER 2201 parameter (MASTER=0) 2202 2203 integer taskid 2204 character*255 full_bak 2205 2206 call Parallel_taskid(taskid) 2207 2208 if (taskid.eq.MASTER) then 2209 close(unit=54) 2210 2211* **** remove backup file *** 2212 call util_file_name_noprefix('TAMD-Z99-bak',.false., 2213 > .false., 2214 > full_bak) 2215 call util_file_unlink(full_bak) 2216 end if 2217 2218 2219 return 2220 end 2221 2222* *********************************************** 2223* * * 2224* * tamd_gauss_write * 2225* * * 2226* *********************************************** 2227 subroutine tamd_gauss_write(n,theta,z,vz,E) 2228 implicit none 2229 integer n 2230 real*8 theta(*),z(*),vz(*),E(*) 2231 2232* **** local variables **** 2233 integer MASTER,taskid 2234 parameter (MASTER=0) 2235 2236 integer i 2237 2238 call Parallel_taskid(taskid) 2239 2240 if (taskid.eq.MASTER) then 2241 write(54,111) (theta(i),i=1,n),(z(i),i=1,n),(vz(i),i=1,n), 2242 > E(2)-E(33)+E(34),E(3),E(4),E(33),E(34) 2243 call util_flush(54) 2244 end if 2245 111 format(99e14.6) 2246 2247 return 2248 end 2249 2250 2251 2252* *********************************************** 2253* * * 2254* * tamd_gauss_analyze * 2255* * * 2256* *********************************************** 2257 subroutine tamd_gauss_analyze(rtdb,filename) 2258 implicit none 2259 integer rtdb 2260 character*(*) filename 2261 2262#include "bafdecls.fh" 2263#include "btdb.fh" 2264#include "stdio.fh" 2265#include "tamd.fh" 2266 2267* **** local variables **** 2268 integer taskid,MASTER 2269 parameter (MASTER=0) 2270 2271 logical found,samez,value 2272 integer i,s1,s2,ndim,nparam,nindx,ntype,count,sindx(10),sparam(10) 2273 integer indx(50) 2274 real*8 gamma,kspring 2275 real*8 f,ff,zz,fold,param(50),tmp(50) 2276 real*8 z(10),k(10),fsum(10),fsum2(10),fdelta(10),zave(10) 2277 character*255 full_filename,atmp 2278 character*500 eqnstring 2279 2280 call Parallel_taskid(taskid) 2281 2282 value = btdb_parallel(.false.) 2283 if ((taskid.eq.MASTER).and.(zfixed)) then 2284 2285 call util_file_name_noprefix(filename,.false., 2286 > .false., 2287 > full_filename) 2288 open(unit=55,file=full_filename,form='formatted', 2289 > status='old') 2290 2291 read(55,*,end=100) 2292 read(55,*,end=100) atmp,ndim 2293 read(55,*,end=100) 2294 read(55,*,end=100) atmp,(sindx(i), i=1,ndim) 2295 read(55,*,end=100) atmp,(sparam(i),i=1,ndim) 2296 read(55,*,end=100) 2297 read(55,*,end=100) 2298 read(55,*,end=100) 2299 read(55,*,end=100) atmp,nindx 2300 read(55,*,end=100) atmp,(indx(i),i=1,nindx) 2301 read(55,*,end=100) atmp,nparam 2302 read(55,*,end=100) atmp,(param(i),i=1,nparam) 2303 read(55,*,end=100) 2304 read(55,*,end=100) 2305 read(55,*,end=100) 2306 read(55,*,end=100) 2307 read(55,*,end=100) 2308 read(55,*,end=100) 2309 2310 do i=1,ndim 2311 fsum(i) = 0.0d0 2312 fsum2(i) = 0.0d0 2313 fdelta(i) = 0.0d0 2314 zave(i) = 0.0d0 2315 k(i) = param(sparam(i)+2) 2316 z(i) = param(sparam(i)+5) 2317 end do 2318 fold = 0.0d0 2319 2320 count = 0 2321 found = .true. 2322 do while (found) 2323 read(55,111,end=100) (tmp(i),i=1,3*ndim+5) 2324 samez = .true. 2325 do i=1,ndim 2326 if (dabs(tmp(ndim+i)-z(i)).gt.1.0d-9) samez=.false. 2327 end do 2328 if (samez) then 2329 do i=1,ndim 2330 if (count.gt.0) fold = fsum(i)/dble(count) 2331 f = -k(i)*(tmp(i)-tmp(ndim+i)) 2332 fsum(i) = fsum(i) + f 2333 fsum2(i) = fsum2(i) + f*f 2334 fdelta(i) = dabs(fsum(i)/dble(count+1) - fold) 2335 zave(i) = zave(i) + tmp(i) 2336 end do 2337 count = count + 1 2338 end if 2339 end do 2340 2341 100 continue 2342 close(unit=55) 2343 2344 write(luout,*) 2345 write(luout,*) 2346 write(luout,*) "================================" 2347 write(luout,*) "==== TAMD Analysis ====" 2348 write(luout,*) "================================" 2349 write(luout,*) 2350 write(luout,*) " == Average Forces ==" 2351 write(luout,*) 2352 2353 write(luout,112) count 2354 write(luout,113) ndim 2355 2356 do i=1,ndim 2357 s1 = sindx(i) +1 2358 s2 = sparam(i)+1 2359 ntype = indx(s1) 2360 2361 f = fsum(i)/dble(count) 2362 ff = fsum2(i)/dble(count) 2363 zz = zave(i)/dble(count) 2364 if (ntype.eq.1) then 2365 write(luout,120) indx(s1+1),indx(s1+2),z(i),zz,k(i), 2366 > f,fdelta(i),ff-f*f 2367 else if (ntype.eq.2) then 2368 write(luout,220) indx(s1+1),indx(s1+2),indx(s1+3), 2369 > z(i),zz,k(i), 2370 > f,fdelta(i),ff-f*f 2371 else if (ntype.eq.3) then 2372 write(luout,320) indx(s1+1),indx(s1+2), 2373 > indx(s1+3),indx(s1+4), 2374 > z(i),zz,k(i), 2375 > f,fdelta(i),ff-f*f 2376 else if (ntype.eq.4) then 2377 write(luout,420) z(i),zz,k(i), 2378 > f,fdelta(i),ff-f*f 2379 else if (ntype.eq.12) then 2380 call nwpw_expression_eqnstring(rtdb,indx(s1+1),eqnstring) 2381 write(luout,720) trim(eqnstring),z(i),zz,k(i), 2382 > f,fdelta(i),ff-f*f 2383 end if 2384 end do 2385 write(luout,*) 2386 write(luout,114) (z(i), i=1,ndim), 2387 > (zave(i)/dble(count), i=1,ndim), 2388 > (fsum(i)/dble(count),i=1,ndim) 2389 write(luout,*) 2390 end if 2391 value = btdb_parallel(.true.) 2392 2393 2394 111 format(99e14.6) 2395 112 format(1x,'frames used =',I8) 2396 113 format(1x,'number of dimensions =',I8) 2397 114 format(1x,'TAMD Force =',99e14.6) 2398 120 format(1x,'atoms=',I5,I5, 2399 > 1x,'dist=',E10.3,' (ave dist=',E10.3,')', 2400 > 1x,'k=',E10.3, 2401 > 3x,'<F>=',E14.6, 2402 > 2x,'Delta <F>=',E14.6, 2403 > 1x,'(<F**2>-<F>**2=',E10.3,')') 2404 220 format(1x,'atoms=',I5,I5,I5, 2405 > 1x,'angle=',E10.3,' (ave angle=',E10.3,')', 2406 > 1x,'k=',E10.3, 2407 > 3x,'<F>=',E14.6, 2408 > 2x,'Delta <F>=',E14.6, 2409 > 1x,'(<F**2>-<F>**2=',E10.3,')') 2410 320 format(1x,'atoms=',I5,I5,I5,I5, 2411 > 1x,'dihedral=',E10.3,' (ave dihedral=',E10.3,')', 2412 > 1x,'k=',E10.3, 2413 > 3x,'<F>=',E14.6, 2414 > 2x,'Delta <F>=',E14.6, 2415 > 1x,'(<F**2>-<F>**2=',E10.3,')') 2416 420 format(1x,'coord_num=',E10.3,' (ave coord_num=',E10.3,')', 2417 > 1x,'k=',E10.3, 2418 > 3x,'<F>=',E14.6, 2419 > 2x,'Delta <F>=',E14.6, 2420 > 1x,'(<F**2>-<F>**2=',E10.3,')') 2421 720 format(1x,'equation=',A, 2422 > 1x,'z=',E10.3,' (zave=',E10.3,')', 2423 > 1x,'k=',E10.3, 2424 > 3x,'<F>=',E14.6, 2425 > 2x,'Delta <F>=',E14.6, 2426 > 1x,'(<F**2>-<F>**2=',E10.3,')') 2427 2428 2429 return 2430 end 2431 2432 2433 2434 2435 2436 2437 2438* *********************************************** 2439* * * 2440* * tamd_setsn2surface * 2441* * * 2442* *********************************************** 2443 subroutine tamd_setsn2surface(sn2) 2444 implicit none 2445 real*8 sn2(*) 2446 2447#include "bafdecls.fh" 2448#include "errquit.fh" 2449#include "tamd.fh" 2450 2451* **** local variables **** 2452 logical sn2found,blaisbunker 2453 integer d1,d2,i,j(10),d 2454 real*8 xp,xr,ww 2455 real*8 er,ep,erp,betar,betap,betarp,xr0,xp0,xrp0 2456 real*8 wr,ar,br,wp,ap,bp,fr,fp,frp,gr,gp 2457 2458 sn2found = dabs(sn2(1)).gt.1.0d-9 2459 if (sn2found) then 2460 2461 blaisbunker = (sn2(1).lt.0.0d0) 2462 2463 if (blaisbunker) then 2464 er = sn2(2) 2465 ep = sn2(3) 2466 erp = sn2(4) 2467 betar = sn2(5) 2468 betap = sn2(6) 2469 betarp = sn2(7) 2470 xr0 = sn2(8) 2471 xp0 = sn2(9) 2472 xrp0 = sn2(9) 2473 ap = sn2(10) 2474 bp = sn2(11) 2475 2476 else 2477 xr0 = sn2(2) 2478 wr = sn2(3) 2479 ar = sn2(4) 2480 br = sn2(5) 2481 2482 xp0 = sn2(6) 2483 wp = sn2(7) 2484 ap = sn2(8) 2485 bp = sn2(9) 2486 2487 er = sn2(10) 2488 ep = sn2(11) 2489 erp = sn2(12) 2490 end if 2491 2492 if (tamdfound) then 2493 if (nmeta.gt.0) then 2494 do d=1,nmeta 2495 j(d) = 0 2496 end do 2497 do i=1,nxmeta_all 2498 ww = 0.0d0 2499 d1 = 1 2500 d2 = 2 2501 xr = dbl_mb(ameta(1)+d1-1) 2502 > + j(d1)*(dbl_mb(bmeta(1)+d1-1)-dbl_mb(ameta(1)+d1-1)) 2503 > /dble(int_mb(nxmeta(1)+d1-1)-1) 2504 xp = dbl_mb(ameta(1)+d2-1) 2505 > + j(d2)*(dbl_mb(bmeta(1)+d2-1)-dbl_mb(ameta(1)+d2-1)) 2506 > /dble(int_mb(nxmeta(1)+d2-1)-1) 2507 2508 if (blaisbunker) then 2509 fr = er - er*(1.0d0-dexp(-betar*(xr-xr0)))**2 2510 fp = ep - ep*(1.0d0-dexp(-betap*(xp-xp0)))**2 2511 > + (1.0d0-dtanh(ap*xr+bp)) 2512 > *(ep - ep*dexp(-betap*(xp-xp0))) 2513 frp = erp - erp*(1.0d0-dexp(-betarp*(xr+xp-xrp0)))**2 2514 ww = fr+fp+frp 2515 else 2516 fr = er + (erp-er)*(1.0d0+dtanh( (xp0-xp)/ar)) 2517 gr = dexp(-((xr-xr0)/wr)**2) 2518 > *(0.5d0+0.5d0*dtanh((xp-xp0)/br)) 2519 2520 fp = ep + (erp-ep)*(1.0d0+dtanh( (xr0-xr)/ap)) 2521 gp = dexp(-((xp-xp0)/wp)**2) 2522 > *(0.5d0+0.5d0*dtanh((xr-xr0)/bp)) 2523 ww = fr*gr + fp*gp 2524 end if 2525 2526 if (ww.gt.0.0d0) ww = 0.0d0 2527 dbl_mb(ymeta(1)+i-1) = dbl_mb(ymeta(1)+i-1) - ww 2528 2529 j(1) = j(1) + 1 2530 do d=1,nmeta-1 2531 if (j(d).ge.int_mb(nxmeta(1)+d-1)) then 2532 j(d) = 0 2533 j(d+1) = j(d+1)+1 2534 end if 2535 end do 2536 end do 2537 end if 2538 end if 2539 end if 2540 return 2541 end 2542 2543