|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00012 MODULE constraint_clv 00013 USE cell_types, ONLY: cell_type 00014 USE colvar_methods, ONLY: colvar_eval_mol_f 00015 USE colvar_types, ONLY: colvar_type,& 00016 diff_colvar 00017 USE f77_blas 00018 USE input_section_types, ONLY: section_vals_get,& 00019 section_vals_get_subs_vals,& 00020 section_vals_type,& 00021 section_vals_val_get,& 00022 section_vals_val_set 00023 USE kinds, ONLY: dp 00024 USE mathlib, ONLY: matvec_3x3 00025 USE molecule_kind_types, ONLY: colvar_constraint_type,& 00026 fixd_constraint_type,& 00027 get_molecule_kind,& 00028 molecule_kind_type 00029 USE molecule_types_new, ONLY: get_molecule,& 00030 global_constraint_type,& 00031 local_colvar_constraint_type,& 00032 molecule_type 00033 USE particle_types, ONLY: particle_type 00034 #include "cp_common_uses.h" 00035 00036 IMPLICIT NONE 00037 00038 PRIVATE 00039 PUBLIC :: shake_roll_colv_int,& 00040 rattle_roll_colv_int,& 00041 shake_colv_int,& 00042 rattle_colv_int,& 00043 shake_roll_colv_ext,& 00044 rattle_roll_colv_ext,& 00045 shake_colv_ext,& 00046 rattle_colv_ext,& 00047 shake_update_colv_int,& 00048 shake_update_colv_ext 00049 00050 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_clv' 00051 00052 CONTAINS 00053 00054 ! ***************************************************************************** 00062 SUBROUTINE shake_colv_int( molecule, particle_set, pos, vel, dt, ishake,& 00063 cell, imass, max_sigma, error ) 00064 00065 TYPE(molecule_type), POINTER :: molecule 00066 TYPE(particle_type), POINTER :: particle_set( : ) 00067 REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : ) 00068 REAL(kind=dp), INTENT(in) :: dt 00069 INTEGER, INTENT(IN) :: ishake 00070 TYPE(cell_type), POINTER :: cell 00071 REAL(KIND=dp), DIMENSION(:) :: imass 00072 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00073 TYPE(cp_error_type), INTENT(inout) :: error 00074 00075 CHARACTER(len=*), PARAMETER :: routineN = 'shake_colv_int', 00076 routineP = moduleN//':'//routineN 00077 00078 LOGICAL :: failure 00079 TYPE(colvar_constraint_type), POINTER :: colv_list( : ) 00080 TYPE(fixd_constraint_type), 00081 DIMENSION(:), POINTER :: fixd_list 00082 TYPE(local_colvar_constraint_type), 00083 POINTER :: lcolv( : ) 00084 TYPE(molecule_kind_type), POINTER :: molecule_kind 00085 00086 NULLIFY(fixd_list) 00087 failure = .FALSE. 00088 molecule_kind => molecule % molecule_kind 00089 CALL get_molecule_kind ( molecule_kind, colv_list=colv_list, fixd_list=fixd_list ) 00090 CALL get_molecule ( molecule, lcolv=lcolv ) 00091 ! Real Shake 00092 CALL shake_colv_low( fixd_list, colv_list, lcolv, & 00093 particle_set, pos, vel, dt, ishake, cell, imass, max_sigma,& 00094 error ) 00095 00096 END SUBROUTINE shake_colv_int 00097 00098 ! ***************************************************************************** 00103 SUBROUTINE shake_update_colv_int( molecule, dt, motion_section, error) 00104 00105 TYPE(molecule_type), POINTER :: molecule 00106 REAL(kind=dp), INTENT(in) :: dt 00107 TYPE(section_vals_type), POINTER :: motion_section 00108 TYPE(cp_error_type), INTENT(inout) :: error 00109 00110 CHARACTER(len=*), PARAMETER :: routineN = 'shake_update_colv_int', 00111 routineP = moduleN//':'//routineN 00112 00113 LOGICAL :: failure 00114 TYPE(colvar_constraint_type), POINTER :: colv_list( : ) 00115 TYPE(molecule_kind_type), POINTER :: molecule_kind 00116 00117 failure = .FALSE. 00118 molecule_kind => molecule % molecule_kind 00119 CALL get_molecule_kind ( molecule_kind, colv_list=colv_list) 00120 ! Real update of the Shake target 00121 CALL shake_update_colv_low( colv_list, dt, motion_section, error) 00122 00123 END SUBROUTINE shake_update_colv_int 00124 00125 ! ***************************************************************************** 00133 SUBROUTINE rattle_colv_int( molecule, particle_set, vel, dt, irattle,& 00134 cell, imass, max_sigma, error ) 00135 00136 TYPE(molecule_type), POINTER :: molecule 00137 TYPE(particle_type), POINTER :: particle_set( : ) 00138 REAL(KIND=dp), INTENT(INOUT) :: vel( :, : ) 00139 REAL(kind=dp), INTENT(in) :: dt 00140 INTEGER, INTENT(IN) :: irattle 00141 TYPE(cell_type), POINTER :: cell 00142 REAL(KIND=dp), DIMENSION(:) :: imass 00143 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00144 TYPE(cp_error_type), INTENT(inout) :: error 00145 00146 TYPE(colvar_constraint_type), POINTER :: colv_list( : ) 00147 TYPE(fixd_constraint_type), 00148 DIMENSION(:), POINTER :: fixd_list 00149 TYPE(local_colvar_constraint_type), 00150 POINTER :: lcolv( : ) 00151 TYPE(molecule_kind_type), POINTER :: molecule_kind 00152 00153 NULLIFY(fixd_list) 00154 molecule_kind => molecule % molecule_kind 00155 CALL get_molecule_kind ( molecule_kind, colv_list = colv_list, fixd_list=fixd_list ) 00156 CALL get_molecule ( molecule, lcolv=lcolv ) 00157 ! Real Rattle 00158 CALL rattle_colv_low ( fixd_list, colv_list, lcolv, & 00159 particle_set, vel, dt, irattle, cell, imass, max_sigma, error ) 00160 00161 END SUBROUTINE rattle_colv_int 00162 00163 ! ***************************************************************************** 00171 SUBROUTINE shake_roll_colv_int( molecule, particle_set, pos, vel, r_shake, v_shake, & 00172 dt, ishake, cell, imass, max_sigma, error ) 00173 00174 TYPE(molecule_type), POINTER :: molecule 00175 TYPE(particle_type), POINTER :: particle_set( : ) 00176 REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : ) 00177 REAL(KIND=dp), DIMENSION(:, :), 00178 INTENT(IN) :: r_shake, v_shake 00179 REAL(kind=dp), INTENT(in) :: dt 00180 INTEGER, INTENT(in) :: ishake 00181 TYPE(cell_type), POINTER :: cell 00182 REAL(KIND=dp), DIMENSION(:) :: imass 00183 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00184 TYPE(cp_error_type), INTENT(inout) :: error 00185 00186 TYPE(colvar_constraint_type), POINTER :: colv_list( : ) 00187 TYPE(fixd_constraint_type), 00188 DIMENSION(:), POINTER :: fixd_list 00189 TYPE(local_colvar_constraint_type), 00190 POINTER :: lcolv( : ) 00191 TYPE(molecule_kind_type), POINTER :: molecule_kind 00192 00193 NULLIFY(fixd_list) 00194 molecule_kind => molecule % molecule_kind 00195 CALL get_molecule_kind ( molecule_kind, colv_list = colv_list, fixd_list=fixd_list ) 00196 CALL get_molecule ( molecule, lcolv=lcolv ) 00197 ! Real Shake 00198 CALL shake_roll_colv_low( fixd_list, colv_list, lcolv, & 00199 particle_set, pos, vel, r_shake, v_shake, dt, ishake, cell,& 00200 imass, max_sigma, error ) 00201 00202 END SUBROUTINE shake_roll_colv_int 00203 00204 ! ***************************************************************************** 00212 SUBROUTINE rattle_roll_colv_int ( molecule, particle_set, vel, r_rattle, & 00213 dt, irattle, veps, cell, imass, max_sigma, error ) 00214 00215 TYPE(molecule_type), POINTER :: molecule 00216 TYPE(particle_type), POINTER :: particle_set( : ) 00217 REAL(KIND=dp), INTENT(INOUT) :: vel( :, : ) 00218 REAL(KIND=dp), INTENT(IN) :: r_rattle( :, : ), dt 00219 INTEGER, INTENT(in) :: irattle 00220 REAL(KIND=dp), INTENT(IN) :: veps( :, : ) 00221 TYPE(cell_type), POINTER :: cell 00222 REAL(KIND=dp), DIMENSION(:) :: imass 00223 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00224 TYPE(cp_error_type), INTENT(inout) :: error 00225 00226 TYPE(colvar_constraint_type), POINTER :: colv_list( : ) 00227 TYPE(fixd_constraint_type), 00228 DIMENSION(:), POINTER :: fixd_list 00229 TYPE(local_colvar_constraint_type), 00230 POINTER :: lcolv( : ) 00231 TYPE(molecule_kind_type), POINTER :: molecule_kind 00232 00233 NULLIFY(fixd_list) 00234 molecule_kind => molecule % molecule_kind 00235 CALL get_molecule_kind ( molecule_kind, colv_list = colv_list, fixd_list=fixd_list ) 00236 CALL get_molecule ( molecule, lcolv=lcolv ) 00237 ! Real Rattle 00238 CALL rattle_roll_colv_low (fixd_list, colv_list, lcolv, & 00239 particle_set, vel, r_rattle, dt, irattle, veps, cell,& 00240 imass, max_sigma, error ) 00241 00242 END SUBROUTINE rattle_roll_colv_int 00243 00244 ! ***************************************************************************** 00252 SUBROUTINE shake_colv_ext( gci, particle_set, pos, vel, dt, ishake,& 00253 cell, imass, max_sigma, error ) 00254 00255 TYPE(global_constraint_type), POINTER :: gci 00256 TYPE(particle_type), POINTER :: particle_set( : ) 00257 REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : ) 00258 REAL(kind=dp), INTENT(in) :: dt 00259 INTEGER, INTENT(IN) :: ishake 00260 TYPE(cell_type), POINTER :: cell 00261 REAL(KIND=dp), DIMENSION(:) :: imass 00262 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00263 TYPE(cp_error_type), INTENT(inout) :: error 00264 00265 CHARACTER(len=*), PARAMETER :: routineN = 'shake_colv_ext', 00266 routineP = moduleN//':'//routineN 00267 00268 TYPE(colvar_constraint_type), POINTER :: colv_list( : ) 00269 TYPE(fixd_constraint_type), 00270 DIMENSION(:), POINTER :: fixd_list 00271 TYPE(local_colvar_constraint_type), 00272 POINTER :: lcolv( : ) 00273 00274 colv_list => gci%colv_list 00275 fixd_list => gci%fixd_list 00276 lcolv => gci%lcolv 00277 ! Real Shake 00278 CALL shake_colv_low( fixd_list, colv_list, lcolv, & 00279 particle_set, pos, vel, dt, ishake, cell, imass, max_sigma,& 00280 error ) 00281 00282 END SUBROUTINE shake_colv_ext 00283 00284 ! ***************************************************************************** 00289 SUBROUTINE shake_update_colv_ext( gci, dt, motion_section, error) 00290 00291 TYPE(global_constraint_type), POINTER :: gci 00292 REAL(kind=dp), INTENT(in) :: dt 00293 TYPE(section_vals_type), POINTER :: motion_section 00294 TYPE(cp_error_type), INTENT(inout) :: error 00295 00296 CHARACTER(len=*), PARAMETER :: routineN = 'shake_update_colv_ext', 00297 routineP = moduleN//':'//routineN 00298 00299 TYPE(colvar_constraint_type), POINTER :: colv_list( : ) 00300 00301 colv_list => gci%colv_list 00302 ! Real update of the Shake target 00303 CALL shake_update_colv_low( colv_list, dt, motion_section, error) 00304 00305 END SUBROUTINE shake_update_colv_ext 00306 00307 ! ***************************************************************************** 00315 SUBROUTINE rattle_colv_ext( gci, particle_set, vel, dt, irattle,& 00316 cell, imass, max_sigma, error ) 00317 00318 TYPE(global_constraint_type), POINTER :: gci 00319 TYPE(particle_type), POINTER :: particle_set( : ) 00320 REAL(KIND=dp), INTENT(INOUT) :: vel( :, : ) 00321 REAL(kind=dp), INTENT(in) :: dt 00322 INTEGER, INTENT(IN) :: irattle 00323 TYPE(cell_type), POINTER :: cell 00324 REAL(KIND=dp), DIMENSION(:) :: imass 00325 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00326 TYPE(cp_error_type), INTENT(inout) :: error 00327 00328 TYPE(colvar_constraint_type), POINTER :: colv_list( : ) 00329 TYPE(fixd_constraint_type), 00330 DIMENSION(:), POINTER :: fixd_list 00331 TYPE(local_colvar_constraint_type), 00332 POINTER :: lcolv( : ) 00333 00334 colv_list => gci%colv_list 00335 fixd_list => gci%fixd_list 00336 lcolv => gci%lcolv 00337 ! Real Rattle 00338 CALL rattle_colv_low ( fixd_list, colv_list, lcolv, & 00339 particle_set, vel, dt, irattle, cell, imass, max_sigma, error ) 00340 00341 END SUBROUTINE rattle_colv_ext 00342 00343 ! ***************************************************************************** 00351 SUBROUTINE shake_roll_colv_ext( gci, particle_set, pos, vel, r_shake, v_shake, & 00352 dt, ishake, cell, imass, max_sigma, error ) 00353 00354 TYPE(global_constraint_type), POINTER :: gci 00355 TYPE(particle_type), POINTER :: particle_set( : ) 00356 REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : ) 00357 REAL(KIND=dp), DIMENSION(:, :), 00358 INTENT(IN) :: r_shake, v_shake 00359 REAL(kind=dp), INTENT(in) :: dt 00360 INTEGER, INTENT(in) :: ishake 00361 TYPE(cell_type), POINTER :: cell 00362 REAL(KIND=dp), DIMENSION(:) :: imass 00363 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00364 TYPE(cp_error_type), INTENT(inout) :: error 00365 00366 TYPE(colvar_constraint_type), POINTER :: colv_list( : ) 00367 TYPE(fixd_constraint_type), 00368 DIMENSION(:), POINTER :: fixd_list 00369 TYPE(local_colvar_constraint_type), 00370 POINTER :: lcolv( : ) 00371 00372 colv_list => gci%colv_list 00373 fixd_list => gci%fixd_list 00374 lcolv => gci%lcolv 00375 ! Real Shake 00376 CALL shake_roll_colv_low( fixd_list, colv_list, lcolv, & 00377 particle_set, pos, vel, r_shake, v_shake, dt, ishake, cell,& 00378 imass, max_sigma, error ) 00379 00380 END SUBROUTINE shake_roll_colv_ext 00381 00382 ! ***************************************************************************** 00390 SUBROUTINE rattle_roll_colv_ext ( gci, particle_set, vel, r_rattle, & 00391 dt, irattle, veps, cell, imass, max_sigma, error ) 00392 00393 TYPE(global_constraint_type), POINTER :: gci 00394 TYPE(particle_type), POINTER :: particle_set( : ) 00395 REAL(KIND=dp), INTENT(INOUT) :: vel( :, : ) 00396 REAL(KIND=dp), INTENT(IN) :: r_rattle( :, : ), dt 00397 INTEGER, INTENT(in) :: irattle 00398 REAL(KIND=dp), INTENT(IN) :: veps( :, : ) 00399 TYPE(cell_type), POINTER :: cell 00400 REAL(KIND=dp), DIMENSION(:) :: imass 00401 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00402 TYPE(cp_error_type), INTENT(inout) :: error 00403 00404 TYPE(colvar_constraint_type), POINTER :: colv_list( : ) 00405 TYPE(fixd_constraint_type), 00406 DIMENSION(:), POINTER :: fixd_list 00407 TYPE(local_colvar_constraint_type), 00408 POINTER :: lcolv( : ) 00409 00410 colv_list => gci%colv_list 00411 fixd_list => gci%fixd_list 00412 lcolv => gci%lcolv 00413 ! Real Rattle 00414 CALL rattle_roll_colv_low (fixd_list, colv_list, lcolv, & 00415 particle_set, vel, r_rattle, dt, irattle, veps, cell,& 00416 imass, max_sigma, error ) 00417 00418 END SUBROUTINE rattle_roll_colv_ext 00419 00420 ! ***************************************************************************** 00428 SUBROUTINE shake_colv_low( fixd_list, colv_list, lcolv, & 00429 particle_set, pos, vel, dt, ishake, cell, imass, max_sigma,& 00430 error ) 00431 TYPE(fixd_constraint_type), 00432 DIMENSION(:), POINTER :: fixd_list 00433 TYPE(colvar_constraint_type), POINTER :: colv_list( : ) 00434 TYPE(local_colvar_constraint_type), 00435 POINTER :: lcolv( : ) 00436 TYPE(particle_type), POINTER :: particle_set( : ) 00437 REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : ) 00438 REAL(kind=dp), INTENT(in) :: dt 00439 INTEGER, INTENT(IN) :: ishake 00440 TYPE(cell_type), POINTER :: cell 00441 REAL(KIND=dp), DIMENSION(:) :: imass 00442 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00443 TYPE(cp_error_type), INTENT(inout) :: error 00444 00445 CHARACTER(len=*), PARAMETER :: routineN = 'shake_colv_low', 00446 routineP = moduleN//':'//routineN 00447 00448 INTEGER :: iconst 00449 LOGICAL :: failure 00450 REAL(KIND=dp) :: del_lam, dtby2, dtsqby2, 00451 fdotf_sum 00452 00453 failure = .FALSE. 00454 dtsqby2 = dt*dt*.5_dp 00455 dtby2 = dt*.5_dp 00456 IF (ishake==1) THEN 00457 DO iconst = 1, SIZE(colv_list) 00458 IF (colv_list(iconst)%restraint%active) CYCLE 00459 ! Update positions 00460 CALL update_con_colv(pos, dtsqby2, lcolv(iconst), & 00461 lambda=lcolv(iconst)%lambda,& 00462 imass=imass, error=error) 00463 ! Update velocities 00464 CALL update_con_colv(vel, dtby2, lcolv(iconst), & 00465 lambda=lcolv(iconst)%lambda,& 00466 imass=imass, error=error) 00467 END DO 00468 ELSE 00469 DO iconst = 1, SIZE(colv_list) 00470 IF (colv_list(iconst)%restraint%active) CYCLE 00471 ! Update colvar 00472 CALL colvar_eval_mol_f( lcolv ( iconst ) % colvar, cell, particles=particle_set,& 00473 pos=pos, fixd_list=fixd_list, error=error) 00474 lcolv ( iconst ) % sigma = diff_colvar(lcolv(iconst)%colvar,& 00475 colv_list(iconst)%expected_value) 00476 fdotf_sum = eval_Jac_colvar(lcolv ( iconst ) % colvar,& 00477 lcolv ( iconst ) % colvar_old, imass=imass, error=error) 00478 del_lam = 2.0_dp*lcolv ( iconst ) % sigma/(dt*dt*fdotf_sum) 00479 lcolv ( iconst ) % lambda = lcolv ( iconst ) % lambda + del_lam 00480 00481 ! Update positions 00482 CALL update_con_colv(pos, dtsqby2, lcolv(iconst), & 00483 lambda=del_lam,& 00484 imass=imass, error=error) 00485 ! Update velocities 00486 CALL update_con_colv(vel, dtby2, lcolv(iconst), & 00487 lambda=del_lam,& 00488 imass=imass, error=error) 00489 END DO 00490 END IF 00491 ! computing the constraint and value of tolerance 00492 DO iconst = 1, SIZE(colv_list) 00493 IF (colv_list(iconst)%restraint%active) CYCLE 00494 CALL colvar_eval_mol_f( lcolv ( iconst ) % colvar, cell, particles=particle_set,& 00495 pos=pos, fixd_list=fixd_list, error=error) 00496 lcolv ( iconst ) % sigma = diff_colvar(lcolv ( iconst ) % colvar,& 00497 colv_list ( iconst ) % expected_value) 00498 max_sigma = MAX(ABS(lcolv ( iconst ) % sigma),max_sigma) 00499 END DO 00500 END SUBROUTINE shake_colv_low 00501 00502 ! ***************************************************************************** 00507 SUBROUTINE shake_update_colv_low(colv_list, dt, motion_section, error ) 00508 TYPE(colvar_constraint_type), POINTER :: colv_list( : ) 00509 REAL(kind=dp), INTENT(in) :: dt 00510 TYPE(section_vals_type), POINTER :: motion_section 00511 TYPE(cp_error_type), INTENT(inout) :: error 00512 00513 CHARACTER(len=*), PARAMETER :: routineN = 'shake_update_colv_low', 00514 routineP = moduleN//':'//routineN 00515 00516 INTEGER :: iconst, irep, n_rep 00517 LOGICAL :: do_update_colvar, explicit, 00518 failure 00519 REAL(KIND=dp) :: clv_target, limit, 00520 new_clv_target, value 00521 TYPE(section_vals_type), POINTER :: collective_sections 00522 00523 failure = .FALSE. 00524 ! Update globally for restart 00525 collective_sections => section_vals_get_subs_vals(motion_section,"CONSTRAINT%COLLECTIVE",error=error) 00526 CALL section_vals_get(collective_sections, n_repetition=n_rep, error=error) 00527 IF (n_rep/=0) THEN 00528 DO irep = 1, n_rep 00529 CALL section_vals_val_get(collective_sections, "TARGET_GROWTH", r_val=value,& 00530 i_rep_section=irep, error=error) 00531 IF ( value /= 0.0_dp ) THEN 00532 CALL section_vals_val_get(collective_sections, "TARGET", r_val=clv_target,& 00533 i_rep_section=irep, error=error) 00534 new_clv_target = clv_target + value * dt 00535 ! Check limits.. 00536 CALL section_vals_val_get(collective_sections, "TARGET_LIMIT", explicit=explicit,& 00537 i_rep_section=irep, error=error) 00538 do_update_colvar = .TRUE. 00539 IF (explicit) THEN 00540 CALL section_vals_val_get(collective_sections, "TARGET_LIMIT", r_val=limit,& 00541 i_rep_section=irep, error=error) 00542 IF (value>0.0_dp) THEN 00543 IF (clv_target==limit) THEN 00544 do_update_colvar = .FALSE. 00545 ELSE IF (new_clv_target>=limit) THEN 00546 new_clv_target = limit 00547 END IF 00548 ELSE 00549 IF (clv_target==limit) THEN 00550 do_update_colvar = .FALSE. 00551 ELSE IF (new_clv_target<=limit) THEN 00552 new_clv_target = limit 00553 END IF 00554 ENDIF 00555 END IF 00556 IF (do_update_colvar) THEN 00557 CALL section_vals_val_set(collective_sections, "TARGET", r_val=new_clv_target,& 00558 i_rep_section=irep, error=error) 00559 END IF 00560 END IF 00561 END DO 00562 END IF 00563 00564 ! Update locally the value to each processor 00565 DO iconst = 1, SIZE(colv_list) 00566 ! Update local to each processor 00567 IF (colv_list(iconst)%expected_value_growth_speed == 0.0_dp) CYCLE 00568 CALL section_vals_val_get(collective_sections, "TARGET", & 00569 r_val=colv_list(iconst)%expected_value,& 00570 i_rep_section=colv_list(iconst)%inp_seq_num, error=error) 00571 END DO 00572 00573 END SUBROUTINE shake_update_colv_low 00574 00575 ! ***************************************************************************** 00583 SUBROUTINE rattle_colv_low ( fixd_list, colv_list, lcolv, & 00584 particle_set, vel, dt, irattle, cell, imass, max_sigma, error ) 00585 00586 TYPE(fixd_constraint_type), 00587 DIMENSION(:), POINTER :: fixd_list 00588 TYPE(colvar_constraint_type), POINTER :: colv_list( : ) 00589 TYPE(local_colvar_constraint_type), 00590 POINTER :: lcolv( : ) 00591 TYPE(particle_type), POINTER :: particle_set( : ) 00592 REAL(KIND=dp), INTENT(INOUT) :: vel( :, : ) 00593 REAL(kind=dp), INTENT(in) :: dt 00594 INTEGER, INTENT(IN) :: irattle 00595 TYPE(cell_type), POINTER :: cell 00596 REAL(KIND=dp), DIMENSION(:) :: imass 00597 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00598 TYPE(cp_error_type), INTENT(inout) :: error 00599 00600 INTEGER :: iconst 00601 REAL(KIND=dp) :: del_lam, dtby2, fdotf_sum 00602 00603 dtby2 = dt*.5_dp 00604 IF (irattle==1) THEN 00605 DO iconst = 1, SIZE(colv_list) 00606 IF (colv_list(iconst)%restraint%active) CYCLE 00607 ! Update colvar_old 00608 CALL colvar_eval_mol_f(lcolv ( iconst ) % colvar_old, cell,& 00609 particles=particle_set, fixd_list=fixd_list, error=error) 00610 ! Update velocities 00611 CALL update_con_colv(vel, dtby2, lcolv(iconst), & 00612 lambda=lcolv(iconst)%lambda,& 00613 imass=imass,error=error) 00614 END DO 00615 ELSE 00616 DO iconst = 1, SIZE(colv_list) 00617 IF (colv_list(iconst)%restraint%active) CYCLE 00618 lcolv ( iconst ) % sigma = rattle_con_eval(lcolv ( iconst ) % colvar_old, vel,error=error) 00619 fdotf_sum = eval_Jac_colvar(lcolv ( iconst ) % colvar_old,& 00620 lcolv ( iconst ) % colvar_old, imass=imass,error=error) 00621 del_lam = 2.0_dp*lcolv ( iconst ) % sigma/(dt*fdotf_sum) 00622 lcolv ( iconst ) % lambda = lcolv ( iconst ) % lambda + del_lam 00623 00624 ! Update velocities 00625 CALL update_con_colv(vel, dtby2, lcolv(iconst), & 00626 lambda=del_lam,& 00627 imass=imass,error=error) 00628 END DO 00629 END IF 00630 00631 DO iconst = 1, SIZE(colv_list) 00632 IF (colv_list(iconst)%restraint%active) CYCLE 00633 lcolv ( iconst ) % sigma = rattle_con_eval(lcolv ( iconst ) % colvar_old, vel,error=error) 00634 max_sigma = MAX(ABS(lcolv ( iconst ) % sigma),max_sigma) 00635 END DO 00636 00637 END SUBROUTINE rattle_colv_low 00638 00639 ! ***************************************************************************** 00647 SUBROUTINE shake_roll_colv_low( fixd_list, colv_list, lcolv, & 00648 particle_set, pos, vel, r_shake, v_shake, dt, ishake, cell,& 00649 imass, max_sigma, error ) 00650 00651 TYPE(fixd_constraint_type), 00652 DIMENSION(:), POINTER :: fixd_list 00653 TYPE(colvar_constraint_type), POINTER :: colv_list( : ) 00654 TYPE(local_colvar_constraint_type), 00655 POINTER :: lcolv( : ) 00656 TYPE(particle_type), POINTER :: particle_set( : ) 00657 REAL(KIND=dp), INTENT(INOUT) :: pos( :, : ), vel( :, : ) 00658 REAL(KIND=dp), DIMENSION(:, :), 00659 INTENT(IN) :: r_shake, v_shake 00660 REAL(kind=dp), INTENT(in) :: dt 00661 INTEGER, INTENT(in) :: ishake 00662 TYPE(cell_type), POINTER :: cell 00663 REAL(KIND=dp), DIMENSION(:) :: imass 00664 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00665 TYPE(cp_error_type), INTENT(inout) :: error 00666 00667 INTEGER :: iconst 00668 REAL(KIND=dp) :: del_lam, dtby2, dtsqby2, 00669 fdotf_sum 00670 00671 dtsqby2 = dt*dt*.5_dp 00672 dtby2 = dt*.5_dp 00673 IF (ishake==1) THEN 00674 DO iconst = 1, SIZE(colv_list) 00675 IF (colv_list(iconst)%restraint%active) CYCLE 00676 ! Update positions 00677 CALL update_con_colv(pos, dtsqby2, lcolv(iconst), & 00678 lambda=lcolv(iconst)%lambda,& 00679 roll=.TRUE.,rmat=r_shake,imass=imass,error=error) 00680 ! Update velocities 00681 CALL update_con_colv(vel, dtby2, lcolv(iconst), & 00682 lambda=lcolv(iconst)%lambda,& 00683 roll=.TRUE.,rmat=v_shake,imass=imass,error=error) 00684 END DO 00685 ELSE 00686 DO iconst = 1, SIZE(colv_list) 00687 IF (colv_list(iconst)%restraint%active) CYCLE 00688 ! Update colvar 00689 CALL colvar_eval_mol_f( lcolv ( iconst ) % colvar, cell, particles=particle_set,& 00690 pos=pos, fixd_list=fixd_list, error=error) 00691 lcolv ( iconst ) % sigma = diff_colvar(lcolv ( iconst ) % colvar,& 00692 colv_list ( iconst ) % expected_value) 00693 fdotf_sum = eval_Jac_colvar(lcolv ( iconst ) % colvar,& 00694 lcolv ( iconst ) % colvar_old, roll=.TRUE., rmat=r_shake,& 00695 imass=imass, error=error) 00696 del_lam = 2.0_dp*lcolv ( iconst ) % sigma/(dt*dt*fdotf_sum) 00697 lcolv ( iconst ) % lambda = lcolv ( iconst ) % lambda + del_lam 00698 00699 ! Update positions 00700 CALL update_con_colv(pos, dtsqby2, lcolv(iconst), & 00701 lambda=del_lam,& 00702 roll=.TRUE., rmat=r_shake, imass=imass,error=error) 00703 ! Update velocities 00704 CALL update_con_colv(vel, dtby2, lcolv(iconst), & 00705 lambda=del_lam,& 00706 roll=.TRUE., rmat=v_shake, imass=imass,error=error) 00707 END DO 00708 END IF 00709 ! computing the constraint and value of tolerance 00710 DO iconst = 1, SIZE(colv_list) 00711 IF (colv_list(iconst)%restraint%active) CYCLE 00712 CALL colvar_eval_mol_f( lcolv ( iconst ) % colvar, cell, particles=particle_set,& 00713 pos=pos, fixd_list=fixd_list, error=error) 00714 lcolv ( iconst ) % sigma = diff_colvar(lcolv ( iconst ) % colvar,& 00715 colv_list ( iconst ) % expected_value) 00716 max_sigma = MAX(ABS(lcolv ( iconst ) % sigma),max_sigma) 00717 END DO 00718 00719 END SUBROUTINE shake_roll_colv_low 00720 00721 ! ***************************************************************************** 00729 SUBROUTINE rattle_roll_colv_low (fixd_list, colv_list, lcolv, & 00730 particle_set, vel, r_rattle, dt, irattle, veps, cell,& 00731 imass, max_sigma, error ) 00732 00733 TYPE(fixd_constraint_type), 00734 DIMENSION(:), POINTER :: fixd_list 00735 TYPE(colvar_constraint_type), POINTER :: colv_list( : ) 00736 TYPE(local_colvar_constraint_type), 00737 POINTER :: lcolv( : ) 00738 TYPE(particle_type), POINTER :: particle_set( : ) 00739 REAL(KIND=dp), INTENT(INOUT) :: vel( :, : ) 00740 REAL(KIND=dp), INTENT(IN) :: r_rattle( :, : ), dt 00741 INTEGER, INTENT(in) :: irattle 00742 REAL(KIND=dp), INTENT(IN) :: veps( :, : ) 00743 TYPE(cell_type), POINTER :: cell 00744 REAL(KIND=dp), DIMENSION(:) :: imass 00745 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 00746 TYPE(cp_error_type), INTENT(inout) :: error 00747 00748 INTEGER :: iconst 00749 REAL(KIND=dp) :: del_lam, dtby2, fdotf_sum 00750 00751 dtby2 = dt*.5_dp 00752 IF (irattle==1) THEN 00753 DO iconst = 1, SIZE ( colv_list ) 00754 IF (colv_list(iconst)%restraint%active) CYCLE 00755 ! Update colvar_old 00756 CALL colvar_eval_mol_f(lcolv ( iconst ) % colvar_old, cell,& 00757 particles=particle_set, fixd_list=fixd_list, error=error) 00758 ! Update velocities 00759 CALL update_con_colv(vel, dtby2, lcolv(iconst), & 00760 lambda=lcolv(iconst)%lambda,& 00761 imass=imass,error=error) 00762 END DO 00763 ELSE 00764 DO iconst = 1, SIZE ( colv_list ) 00765 IF (colv_list(iconst)%restraint%active) CYCLE 00766 lcolv ( iconst ) % sigma = rattle_con_eval(lcolv ( iconst ) % colvar_old, vel,& 00767 roll=.TRUE., veps=veps, rmat=r_rattle, particles=particle_set,& 00768 error=error) 00769 fdotf_sum = eval_Jac_colvar(lcolv ( iconst ) % colvar_old,& 00770 lcolv ( iconst ) % colvar_old, roll=.TRUE.,& 00771 rmat=r_rattle,imass=imass, error=error) 00772 del_lam = 2.0_dp*lcolv ( iconst ) % sigma/(dt*fdotf_sum) 00773 lcolv ( iconst ) % lambda = lcolv ( iconst ) % lambda + del_lam 00774 ! Update velocities 00775 CALL update_con_colv(vel, dtby2, lcolv(iconst), & 00776 lambda=del_lam,& 00777 roll=.TRUE.,rmat=r_rattle,imass=imass,error=error) 00778 END DO 00779 END IF 00780 ! computing the constraint and value of the tolerance 00781 DO iconst = 1, SIZE(colv_list) 00782 IF (colv_list(iconst)%restraint%active) CYCLE 00783 lcolv ( iconst ) % sigma = rattle_con_eval(lcolv ( iconst ) % colvar_old, vel,& 00784 roll=.TRUE., veps=veps, rmat=r_rattle, particles=particle_set,& 00785 error=error) 00786 max_sigma = MAX(ABS(lcolv ( iconst ) % sigma),max_sigma) 00787 END DO 00788 00789 END SUBROUTINE rattle_roll_colv_low 00790 00791 ! ***************************************************************************** 00797 SUBROUTINE update_con_colv(wrk, fac, lcolv, lambda, roll, rmat, imass, error) 00798 REAL(KIND=dp), INTENT(INOUT) :: wrk( :, : ) 00799 REAL(KIND=dp), INTENT(IN) :: fac 00800 TYPE(local_colvar_constraint_type), 00801 INTENT(IN) :: lcolv 00802 REAL(KIND=dp), INTENT(IN) :: lambda 00803 LOGICAL, INTENT(in), OPTIONAL :: roll 00804 REAL(KIND=dp), DIMENSION(:, :), 00805 INTENT(IN), OPTIONAL :: rmat 00806 REAL(KIND=dp), DIMENSION(:) :: imass 00807 TYPE(cp_error_type), INTENT(inout) :: error 00808 00809 CHARACTER(len=*), PARAMETER :: routineN = 'update_con_colv', 00810 routineP = moduleN//':'//routineN 00811 00812 INTEGER :: iatm, ind 00813 LOGICAL :: failure, my_roll 00814 REAL(KIND=dp), DIMENSION(3) :: f_roll 00815 00816 failure = .FALSE. 00817 my_roll = .FALSE. 00818 IF (PRESENT(roll)) THEN 00819 my_roll = roll 00820 IF (my_roll) THEN 00821 CPPostcondition(PRESENT(rmat),cp_failure_level,routineP,error,failure) 00822 END IF 00823 END IF 00824 IF (.NOT.failure) THEN 00825 DO iatm = 1, SIZE( lcolv % colvar_old % i_atom ) 00826 ind = lcolv % colvar_old % i_atom (iatm) 00827 ! 00828 IF (my_roll) THEN 00829 ! If ROLL rotate forces 00830 CALL matvec_3x3 (f_roll, rmat, lcolv % colvar_old % dsdr(:,iatm) ) 00831 ELSE 00832 f_roll = lcolv % colvar_old % dsdr(:,iatm) 00833 END IF 00834 wrk ( :, ind ) = wrk ( :, ind ) - imass(ind) * fac * lambda * f_roll 00835 END DO 00836 END IF 00837 END SUBROUTINE update_con_colv 00838 00839 ! ***************************************************************************** 00845 FUNCTION eval_Jac_colvar(colvar, colvar_old, roll, rmat, imass, error) RESULT(res) 00846 TYPE(colvar_type), POINTER :: colvar, colvar_old 00847 LOGICAL, INTENT(IN), OPTIONAL :: roll 00848 REAL(KIND=dp), DIMENSION(:, :), 00849 INTENT(IN), OPTIONAL :: rmat 00850 REAL(KIND=dp), DIMENSION(:) :: imass 00851 TYPE(cp_error_type), INTENT(inout) :: error 00852 REAL(KIND=dp) :: res 00853 00854 CHARACTER(len=*), PARAMETER :: routineN = 'eval_Jac_colvar', 00855 routineP = moduleN//':'//routineN 00856 00857 INTEGER :: i, iatom 00858 LOGICAL :: failure, my_roll 00859 REAL(KIND=dp), DIMENSION(3) :: tmp1, tmp2, tmp3 00860 00861 failure = .FALSE. 00862 my_roll = .FALSE. 00863 IF (PRESENT(roll)) THEN 00864 my_roll = roll 00865 IF (my_roll) THEN 00866 CPPostcondition(PRESENT(rmat),cp_failure_level,routineP,error,failure) 00867 END IF 00868 END IF 00869 00870 res = 0.0_dp 00871 DO i = 1, SIZE(colvar%i_atom) 00872 iatom = colvar%i_atom(i) 00873 tmp1 = colvar%dsdr(1:3,i) 00874 IF (my_roll) THEN 00875 tmp3 = colvar_old%dsdr(1:3,i) 00876 CALL matvec_3x3 (tmp2, rmat, tmp3 ) 00877 ELSE 00878 tmp2 = colvar_old%dsdr(1:3,i) 00879 END IF 00880 res = res + DOT_PRODUCT(tmp1,tmp2)*imass(iatom) 00881 END DO 00882 00883 END FUNCTION eval_Jac_colvar 00884 00885 ! ***************************************************************************** 00891 FUNCTION rattle_con_eval(colvar, vel, roll, veps, rmat, particles, error) RESULT(res) 00892 TYPE(colvar_type), POINTER :: colvar 00893 REAL(KIND=dp), INTENT(INOUT) :: vel( :, : ) 00894 LOGICAL, INTENT(IN), OPTIONAL :: roll 00895 REAL(KIND=dp), DIMENSION(:, :), 00896 INTENT(IN), OPTIONAL :: veps, rmat 00897 TYPE(particle_type), DIMENSION(:), 00898 OPTIONAL, POINTER :: particles 00899 TYPE(cp_error_type), INTENT(inout) :: error 00900 REAL(KIND=dp) :: res 00901 00902 CHARACTER(len=*), PARAMETER :: routineN = 'rattle_con_eval', 00903 routineP = moduleN//':'//routineN 00904 00905 INTEGER :: iatm, ind 00906 LOGICAL :: failure, my_roll 00907 REAL(KIND=dp) :: tmp(3) 00908 REAL(KIND=dp), DIMENSION(3) :: f_roll, pos, v_roll 00909 00910 failure = .FALSE. 00911 my_roll = .FALSE. 00912 IF (PRESENT(roll)) THEN 00913 my_roll = roll 00914 IF (my_roll) THEN 00915 CPPostcondition(PRESENT(rmat),cp_failure_level,routineP,error,failure) 00916 CPPostcondition(PRESENT(veps),cp_failure_level,routineP,error,failure) 00917 CPPostcondition(PRESENT(particles),cp_failure_level,routineP,error,failure) 00918 END IF 00919 END IF 00920 res = 0.0_dp 00921 DO iatm = 1, SIZE( colvar % i_atom ) 00922 ind = colvar % i_atom (iatm) 00923 IF (my_roll) THEN 00924 pos = particles ( ind ) % r 00925 CALL matvec_3x3 (f_roll, rmat, colvar % dsdr(:,iatm) ) 00926 CALL matvec_3x3 (tmp, veps, pos ) 00927 v_roll = vel(:,ind) + tmp 00928 ELSE 00929 f_roll = colvar % dsdr(:,iatm) 00930 v_roll = vel(:,ind) 00931 END IF 00932 res = res + DOT_PRODUCT(f_roll, v_roll) 00933 END DO 00934 00935 END FUNCTION rattle_con_eval 00936 00937 END MODULE constraint_clv
1.7.3