|
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 ! ***************************************************************************** 00010 MODULE force_env_utils 00011 00012 USE atomic_kind_list_types, ONLY: atomic_kind_list_type 00013 USE cell_types, ONLY: cell_type 00014 USE constraint, ONLY: rattle_control,& 00015 shake_control 00016 USE constraint_util, ONLY: getold 00017 USE cp_subsys_types, ONLY: cp_subsys_get,& 00018 cp_subsys_type 00019 USE distribution_1d_types, ONLY: distribution_1d_type 00020 USE f77_blas 00021 USE force_env_types, ONLY: force_env_get,& 00022 force_env_type 00023 USE input_section_types, ONLY: section_vals_get,& 00024 section_vals_get_subs_vals,& 00025 section_vals_type,& 00026 section_vals_val_get 00027 USE kinds, ONLY: dp 00028 USE mathlib, ONLY: det_3x3,& 00029 jacobi 00030 USE mol_kind_new_list_types, ONLY: mol_kind_new_list_type 00031 USE mol_new_list_types, ONLY: mol_new_list_type 00032 USE molecule_types_new, ONLY: global_constraint_type 00033 USE particle_list_types, ONLY: particle_list_type 00034 USE particle_types, ONLY: update_particle_set 00035 USE physcon, ONLY: pascal 00036 USE timings, ONLY: timeset,& 00037 timestop 00038 #include "cp_common_uses.h" 00039 00040 IMPLICIT NONE 00041 00042 PRIVATE 00043 00044 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_env_utils' 00045 00046 PUBLIC :: force_env_shake,& 00047 force_env_rattle,& 00048 rescale_forces,& 00049 write_stress_tensor,& 00050 write_forces 00051 00052 CONTAINS 00053 00054 ! ***************************************************************************** 00066 SUBROUTINE force_env_shake(force_env,dt,shake_tol,log_unit,lagrange_mult,dump_lm,& 00067 pos,vel,compold,reset,error) 00068 00069 TYPE(force_env_type), POINTER :: force_env 00070 REAL(kind=dp), INTENT(IN), OPTIONAL :: dt 00071 REAL(kind=dp), INTENT(IN) :: shake_tol 00072 INTEGER, INTENT(in), OPTIONAL :: log_unit, lagrange_mult 00073 LOGICAL, INTENT(IN), OPTIONAL :: dump_lm 00074 REAL(KIND=dp), DIMENSION(:, :), 00075 INTENT(INOUT), OPTIONAL, TARGET :: pos, vel 00076 LOGICAL, INTENT(IN), OPTIONAL :: compold, reset 00077 TYPE(cp_error_type), INTENT(inout) :: error 00078 00079 CHARACTER(len=*), PARAMETER :: routineN = 'force_env_shake', 00080 routineP = moduleN//':'//routineN 00081 00082 INTEGER :: handle, i, iparticle, iparticle_kind, iparticle_local, j, 00083 my_lagrange_mult, my_log_unit, nparticle_kind, nparticle_local, stat 00084 LOGICAL :: failure, has_pos, has_vel, 00085 my_dump_lm 00086 REAL(KIND=dp) :: mydt 00087 REAL(KIND=dp), DIMENSION(:, :), POINTER :: my_pos, my_vel 00088 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 00089 TYPE(cell_type), POINTER :: cell 00090 TYPE(cp_subsys_type), POINTER :: subsys 00091 TYPE(distribution_1d_type), POINTER :: local_molecules, 00092 local_particles 00093 TYPE(global_constraint_type), POINTER :: gci 00094 TYPE(mol_kind_new_list_type), POINTER :: molecule_kinds 00095 TYPE(mol_new_list_type), POINTER :: molecules 00096 TYPE(particle_list_type), POINTER :: particles 00097 00098 failure=.FALSE. 00099 CALL timeset(routineN,handle) 00100 CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) 00101 CPPrecondition(force_env%ref_count>0,cp_failure_level,routineP,error,failure) 00102 IF (.NOT.failure) THEN 00103 my_log_unit=-1 00104 IF (PRESENT(log_unit)) my_log_unit=log_unit 00105 my_lagrange_mult=-1 00106 IF (PRESENT(lagrange_mult)) my_lagrange_mult=lagrange_mult 00107 my_dump_lm = .FALSE. 00108 IF (PRESENT(dump_lm)) my_dump_lm = dump_lm 00109 NULLIFY(subsys,cell,molecules,molecule_kinds,local_molecules,particles,& 00110 my_pos,my_vel,gci) 00111 IF (PRESENT(pos)) my_pos => pos 00112 IF (PRESENT(vel)) my_vel => vel 00113 mydt = 0.1_dp 00114 IF (PRESENT(dt)) mydt = dt 00115 CALL force_env_get(force_env,subsys=subsys,cell=cell,error=error) 00116 CALL cp_subsys_get(subsys, & 00117 atomic_kinds=atomic_kinds,& 00118 local_molecules_new=local_molecules,& 00119 local_particles=local_particles,& 00120 molecules_new=molecules,& 00121 molecule_kinds_new=molecule_kinds,& 00122 particles=particles,& 00123 gci=gci,& 00124 error=error) 00125 nparticle_kind = atomic_kinds%n_els 00126 IF (PRESENT(compold)) THEN 00127 IF (compold) THEN 00128 CALL getold( gci, local_molecules, molecules%els, molecule_kinds%els,& 00129 particles%els, cell, error) 00130 END IF 00131 END IF 00132 has_pos=.FALSE. 00133 IF (.NOT.ASSOCIATED(my_pos)) THEN 00134 has_pos=.TRUE. 00135 ALLOCATE(my_pos(3,particles%n_els),stat=stat) 00136 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00137 my_pos = 0.0_dp 00138 DO iparticle_kind=1,nparticle_kind 00139 nparticle_local = local_particles%n_el(iparticle_kind) 00140 DO iparticle_local=1,nparticle_local 00141 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 00142 my_pos (:,iparticle) = particles%els(iparticle)%r(:) 00143 END DO 00144 END DO 00145 END IF 00146 has_vel=.FALSE. 00147 IF (.NOT.ASSOCIATED(my_vel)) THEN 00148 has_vel=.TRUE. 00149 ALLOCATE(my_vel(3,particles%n_els),stat=stat) 00150 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00151 my_vel = 0.0_dp 00152 DO iparticle_kind=1,nparticle_kind 00153 nparticle_local = local_particles%n_el(iparticle_kind) 00154 DO iparticle_local=1,nparticle_local 00155 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 00156 my_vel (:,iparticle) = particles%els(iparticle)%v(:) 00157 END DO 00158 END DO 00159 END IF 00160 00161 CALL shake_control( gci=gci, local_molecules=local_molecules,& 00162 molecule_set=molecules%els, molecule_kind_set=molecule_kinds%els, & 00163 particle_set=particles%els, pos=my_pos, vel=my_vel, dt=mydt,& 00164 shake_tol=shake_tol, log_unit=my_log_unit, lagrange_mult=my_lagrange_mult,& 00165 dump_lm= my_dump_lm, cell=cell,group=force_env%para_env%group,& 00166 local_particles=local_particles, error=error ) 00167 00168 ! Possibly reset the lagrange multipliers 00169 IF (PRESENT(reset)) THEN 00170 IF (reset) THEN 00171 ! Reset Intramolecular constraints 00172 DO i = 1,SIZE(molecules%els) 00173 IF (ASSOCIATED(molecules%els(i)%lci%lcolv)) THEN 00174 DO j = 1, SIZE(molecules%els(i)%lci%lcolv) 00175 ! Reset langrange multiplier 00176 molecules%els(i)%lci%lcolv(j)%lambda = 0.0_dp 00177 END DO 00178 END IF 00179 IF (ASSOCIATED(molecules%els(i)%lci%lg3x3)) THEN 00180 DO j = 1, SIZE(molecules%els(i)%lci%lg3x3) 00181 ! Reset langrange multiplier 00182 molecules%els(i)%lci%lg3x3(j)%lambda = 0.0_dp 00183 END DO 00184 END IF 00185 IF (ASSOCIATED(molecules%els(i)%lci%lg4x6)) THEN 00186 DO j = 1, SIZE(molecules%els(i)%lci%lg4x6) 00187 ! Reset langrange multiplier 00188 molecules%els(i)%lci%lg4x6(j)%lambda = 0.0_dp 00189 END DO 00190 END IF 00191 END DO 00192 ! Reset Intermolecular constraints 00193 IF (ASSOCIATED(gci)) THEN 00194 IF (ASSOCIATED(gci%lcolv)) THEN 00195 DO j = 1, SIZE(gci%lcolv) 00196 ! Reset langrange multiplier 00197 gci%lcolv(j)%lambda = 0.0_dp 00198 END DO 00199 END IF 00200 IF (ASSOCIATED(gci%lg3x3)) THEN 00201 DO j = 1, SIZE(gci%lg3x3) 00202 ! Reset langrange multiplier 00203 gci%lg3x3(j)%lambda = 0.0_dp 00204 END DO 00205 END IF 00206 IF (ASSOCIATED(gci%lg4x6)) THEN 00207 DO j = 1, SIZE(gci%lg4x6) 00208 ! Reset langrange multiplier 00209 gci%lg4x6(j)%lambda = 0.0_dp 00210 END DO 00211 END IF 00212 END IF 00213 END IF 00214 END IF 00215 00216 IF (has_pos) THEN 00217 CALL update_particle_set ( particles%els, force_env%para_env%group, pos=my_pos,& 00218 error=error) 00219 DEALLOCATE(my_pos,stat=stat) 00220 CPPostcondition(stat==0,cp_warning_level,routineP,error,failure) 00221 END IF 00222 IF (has_vel) THEN 00223 CALL update_particle_set ( particles%els, force_env%para_env%group, vel=my_vel,& 00224 error=error) 00225 DEALLOCATE(my_vel,stat=stat) 00226 CPPostcondition(stat==0,cp_warning_level,routineP,error,failure) 00227 END IF 00228 END IF 00229 CALL timestop(handle) 00230 END SUBROUTINE force_env_shake 00231 00232 ! ***************************************************************************** 00246 SUBROUTINE force_env_rattle(force_env,dt,shake_tol,log_unit,lagrange_mult,dump_lm,& 00247 vel,reset,error) 00248 00249 TYPE(force_env_type), POINTER :: force_env 00250 REAL(kind=dp), INTENT(in), OPTIONAL :: dt 00251 REAL(kind=dp), INTENT(in) :: shake_tol 00252 INTEGER, INTENT(in), OPTIONAL :: log_unit, lagrange_mult 00253 LOGICAL, INTENT(IN), OPTIONAL :: dump_lm 00254 REAL(KIND=dp), DIMENSION(:, :), 00255 INTENT(INOUT), OPTIONAL, TARGET :: vel 00256 LOGICAL, INTENT(IN), OPTIONAL :: reset 00257 TYPE(cp_error_type), INTENT(inout) :: error 00258 00259 CHARACTER(len=*), PARAMETER :: routineN = 'force_env_rattle', 00260 routineP = moduleN//':'//routineN 00261 00262 INTEGER :: handle, i, iparticle, iparticle_kind, iparticle_local, j, 00263 my_lagrange_mult, my_log_unit, nparticle_kind, nparticle_local, stat 00264 LOGICAL :: failure, has_vel, my_dump_lm 00265 REAL(KIND=dp) :: mydt 00266 REAL(KIND=dp), DIMENSION(:, :), POINTER :: my_vel 00267 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 00268 TYPE(cell_type), POINTER :: cell 00269 TYPE(cp_subsys_type), POINTER :: subsys 00270 TYPE(distribution_1d_type), POINTER :: local_molecules, 00271 local_particles 00272 TYPE(global_constraint_type), POINTER :: gci 00273 TYPE(mol_kind_new_list_type), POINTER :: molecule_kinds 00274 TYPE(mol_new_list_type), POINTER :: molecules 00275 TYPE(particle_list_type), POINTER :: particles 00276 00277 failure=.FALSE. 00278 00279 CALL timeset(routineN,handle) 00280 CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) 00281 CPPrecondition(force_env%ref_count>0,cp_failure_level,routineP,error,failure) 00282 IF (.NOT.failure) THEN 00283 my_log_unit=-1 00284 IF (PRESENT(log_unit)) my_log_unit=log_unit 00285 my_lagrange_mult=-1 00286 IF (PRESENT(lagrange_mult)) my_lagrange_mult=lagrange_mult 00287 my_dump_lm = .FALSE. 00288 IF (PRESENT(dump_lm)) my_dump_lm = dump_lm 00289 NULLIFY(subsys,cell,molecules,molecule_kinds,local_molecules,particles,& 00290 my_vel) 00291 IF (PRESENT(vel)) my_vel => vel 00292 mydt = 0.1_dp 00293 IF (PRESENT(dt)) mydt = dt 00294 CALL force_env_get(force_env,subsys=subsys,cell=cell,error=error) 00295 CALL cp_subsys_get(subsys, & 00296 atomic_kinds=atomic_kinds,& 00297 local_molecules_new=local_molecules,& 00298 local_particles=local_particles,& 00299 molecules_new=molecules,& 00300 molecule_kinds_new=molecule_kinds,& 00301 particles=particles,& 00302 gci=gci,& 00303 error=error) 00304 nparticle_kind = atomic_kinds%n_els 00305 has_vel=.FALSE. 00306 IF (.NOT.ASSOCIATED(my_vel)) THEN 00307 has_vel=.TRUE. 00308 ALLOCATE(my_vel(3,particles%n_els),stat=stat) 00309 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00310 my_vel = 0.0_dp 00311 DO iparticle_kind=1,nparticle_kind 00312 nparticle_local = local_particles%n_el(iparticle_kind) 00313 DO iparticle_local=1,nparticle_local 00314 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 00315 my_vel (:,iparticle) = particles%els(iparticle)%v(:) 00316 END DO 00317 END DO 00318 END IF 00319 00320 CALL rattle_control( gci=gci, local_molecules=local_molecules,& 00321 molecule_set=molecules%els, molecule_kind_set=molecule_kinds%els, & 00322 particle_set=particles%els, vel=my_vel, dt=mydt,& 00323 rattle_tol=shake_tol, log_unit=my_log_unit, lagrange_mult=my_lagrange_mult,& 00324 dump_lm=my_dump_lm, cell=cell, group=force_env%para_env%group,& 00325 local_particles=local_particles, error=error ) 00326 00327 ! Possibly reset the lagrange multipliers 00328 IF (PRESENT(reset)) THEN 00329 IF (reset) THEN 00330 ! Reset Intramolecular constraints 00331 DO i = 1,SIZE(molecules%els) 00332 IF (ASSOCIATED(molecules%els(i)%lci%lcolv)) THEN 00333 DO j = 1, SIZE(molecules%els(i)%lci%lcolv) 00334 ! Reset langrange multiplier 00335 molecules%els(i)%lci%lcolv(j)%lambda = 0.0_dp 00336 END DO 00337 END IF 00338 IF (ASSOCIATED(molecules%els(i)%lci%lg3x3)) THEN 00339 DO j = 1, SIZE(molecules%els(i)%lci%lg3x3) 00340 ! Reset langrange multiplier 00341 molecules%els(i)%lci%lg3x3(j)%lambda = 0.0_dp 00342 END DO 00343 END IF 00344 IF (ASSOCIATED(molecules%els(i)%lci%lg4x6)) THEN 00345 DO j = 1, SIZE(molecules%els(i)%lci%lg4x6) 00346 ! Reset langrange multiplier 00347 molecules%els(i)%lci%lg4x6(j)%lambda = 0.0_dp 00348 END DO 00349 END IF 00350 END DO 00351 ! Reset Intermolecular constraints 00352 IF (ASSOCIATED(gci)) THEN 00353 IF (ASSOCIATED(gci%lcolv)) THEN 00354 DO j = 1, SIZE(gci%lcolv) 00355 ! Reset langrange multiplier 00356 gci%lcolv(j)%lambda = 0.0_dp 00357 END DO 00358 END IF 00359 IF (ASSOCIATED(gci%lg3x3)) THEN 00360 DO j = 1, SIZE(gci%lg3x3) 00361 ! Reset langrange multiplier 00362 gci%lg3x3(j)%lambda = 0.0_dp 00363 END DO 00364 END IF 00365 IF (ASSOCIATED(gci%lg4x6)) THEN 00366 DO j = 1, SIZE(gci%lg4x6) 00367 ! Reset langrange multiplier 00368 gci%lg4x6(j)%lambda = 0.0_dp 00369 END DO 00370 END IF 00371 END IF 00372 END IF 00373 END IF 00374 00375 IF (has_vel) THEN 00376 CALL update_particle_set ( particles%els, force_env%para_env%group, vel=my_vel,& 00377 error=error) 00378 END IF 00379 DEALLOCATE(my_vel,stat=stat) 00380 CPPostcondition(stat==0,cp_warning_level,routineP,error,failure) 00381 END IF 00382 CALL timestop(handle) 00383 END SUBROUTINE force_env_rattle 00384 00385 ! ***************************************************************************** 00392 SUBROUTINE rescale_forces (force_env, error) 00393 TYPE(force_env_type), POINTER :: force_env 00394 TYPE(cp_error_type), INTENT(inout) :: error 00395 00396 CHARACTER(len=*), PARAMETER :: routineN = 'rescale_forces', 00397 routineP = moduleN//':'//routineN 00398 00399 INTEGER :: handle, iparticle 00400 LOGICAL :: explicit, failure 00401 REAL(KIND=dp) :: force(3), max_value, mod_force 00402 TYPE(cp_subsys_type), POINTER :: subsys 00403 TYPE(particle_list_type), POINTER :: particles 00404 TYPE(section_vals_type), POINTER :: rescale_force_section 00405 00406 failure=.FALSE. 00407 CALL timeset(routineN,handle) 00408 CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) 00409 CPPrecondition(force_env%ref_count>0,cp_failure_level,routineP,error,failure) 00410 rescale_force_section => section_vals_get_subs_vals(force_env%force_env_section,"RESCALE_FORCES",error=error) 00411 CALL section_vals_get(rescale_force_section, explicit=explicit, error=error) 00412 IF (.NOT.failure.AND.explicit) THEN 00413 CALL section_vals_val_get(rescale_force_section,"MAX_FORCE",r_val=max_value,error=error) 00414 CALL force_env_get(force_env,subsys=subsys,error=error) 00415 CALL cp_subsys_get(subsys,particles=particles,error=error) 00416 DO iparticle = 1, SIZE(particles%els) 00417 force = particles%els(iparticle)%f(:) 00418 mod_force = SQRT(DOT_PRODUCT(force,force)) 00419 IF ((mod_force > max_value).AND.(mod_force /= 0.0_dp)) THEN 00420 force = force / mod_force * max_value 00421 particles%els(iparticle)%f(:) = force 00422 END IF 00423 END DO 00424 END IF 00425 CALL timestop(handle) 00426 END SUBROUTINE rescale_forces 00427 00428 ! ***************************************************************************** 00433 SUBROUTINE write_stress_tensor(pv_virial,output_unit,cell,ndigits,numerical,& 00434 error) 00435 00436 REAL(KIND=dp), DIMENSION(3, 3), 00437 INTENT(IN) :: pv_virial 00438 INTEGER, INTENT(IN) :: output_unit 00439 TYPE(cell_type), POINTER :: cell 00440 INTEGER, INTENT(IN) :: ndigits 00441 LOGICAL, INTENT(IN) :: numerical 00442 TYPE(cp_error_type), INTENT(INOUT) :: error 00443 00444 CHARACTER(LEN=*), PARAMETER :: routineN = 'write_stress_tensor', 00445 routineP = moduleN//':'//routineN 00446 00447 CHARACTER(LEN=15) :: fmtstr3 00448 CHARACTER(LEN=16) :: fmtstr4 00449 CHARACTER(LEN=22) :: fmtstr2 00450 CHARACTER(LEN=27) :: fmtstr5 00451 CHARACTER(LEN=31) :: fmtstr1 00452 INTEGER :: n 00453 LOGICAL :: failure 00454 REAL(KIND=dp), DIMENSION(3) :: eigval 00455 REAL(KIND=dp), DIMENSION(3, 3) :: eigvec, stress_tensor 00456 00457 failure = .FALSE. 00458 00459 IF (output_unit > 0) THEN 00460 CPPrecondition(ASSOCIATED(cell),cp_failure_level,routineP,error,failure) 00461 stress_tensor(:,:) = pv_virial(:,:)/cell%deth*pascal*1.0E-9_dp 00462 n = MIN(MAX(1,ndigits),20) 00463 fmtstr1 = "(/,T2,A,/,/,T13,A1,2( X,A1))" 00464 WRITE (UNIT=fmtstr1(22:23),FMT="(I2)") n + 7 00465 fmtstr2 = "(T3,A,T5,3(1X,F . ))" 00466 WRITE (UNIT=fmtstr2(16:17),FMT="(I2)") n + 7 00467 WRITE (UNIT=fmtstr2(19:20),FMT="(I2)") n 00468 fmtstr3 = "(/,T3,A,F . )" 00469 WRITE (UNIT=fmtstr3(10:11),FMT="(I2)") n + 8 00470 WRITE (UNIT=fmtstr3(13:14),FMT="(I2)") n 00471 IF (numerical) THEN 00472 WRITE (UNIT=output_unit,FMT=fmtstr1)& 00473 "NUMERICAL STRESS TENSOR [GPa]","X","Y","Z" 00474 ELSE 00475 WRITE (UNIT=output_unit,FMT=fmtstr1)& 00476 "STRESS TENSOR [GPa]","X","Y","Z" 00477 END IF 00478 WRITE (UNIT=output_unit,FMT=fmtstr2) "X",stress_tensor(1,1:3) 00479 WRITE (UNIT=output_unit,FMT=fmtstr2) "Y",stress_tensor(2,1:3) 00480 WRITE (UNIT=output_unit,FMT=fmtstr2) "Z",stress_tensor(3,1:3) 00481 fmtstr4 = "(/,T3,A,ES . )" 00482 WRITE (UNIT=fmtstr4(11:12),FMT="(I2)") n + 8 00483 WRITE (UNIT=fmtstr4(14:15),FMT="(I2)") n 00484 WRITE (UNIT=output_unit,FMT=fmtstr4)& 00485 "1/3 Trace(stress tensor): ",(stress_tensor(1,1) +& 00486 stress_tensor(2,2) +& 00487 stress_tensor(3,3))/3.0_dp,& 00488 "Det(stress tensor) : ",det_3x3(stress_tensor(:,1),& 00489 stress_tensor(:,2),& 00490 stress_tensor(:,3)) 00491 eigval(:) = 0.0_dp 00492 eigvec(:,:) = 0.0_dp 00493 CALL jacobi(stress_tensor,eigval,eigvec) 00494 fmtstr5 = "(/,/,T2,A,/,/,T5,3F . ,/)" 00495 WRITE (UNIT=fmtstr5(20:21),FMT="(I2)") n + 8 00496 WRITE (UNIT=fmtstr5(23:24),FMT="(I2)") n 00497 WRITE (UNIT=output_unit,FMT=fmtstr5)& 00498 "EIGENVECTORS AND EIGENVALUES OF THE STRESS TENSOR",& 00499 eigval(1:3) 00500 WRITE (UNIT=output_unit,FMT=fmtstr2) " ",eigvec(1,1:3) 00501 WRITE (UNIT=output_unit,FMT=fmtstr2) " ",eigvec(2,1:3) 00502 WRITE (UNIT=output_unit,FMT=fmtstr2) " ",eigvec(3,1:3) 00503 END IF 00504 00505 END SUBROUTINE write_stress_tensor 00506 00507 ! ***************************************************************************** 00512 SUBROUTINE write_forces(particles,output_unit,label,ndigits,total_force,& 00513 grand_total_force,error) 00514 00515 TYPE(particle_list_type), POINTER :: particles 00516 INTEGER, INTENT(IN) :: output_unit 00517 CHARACTER(LEN=*), INTENT(IN) :: label 00518 INTEGER, INTENT(IN) :: ndigits 00519 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: total_force 00520 REAL(KIND=dp), DIMENSION(3), 00521 INTENT(INOUT), OPTIONAL :: grand_total_force 00522 TYPE(cp_error_type), INTENT(INOUT) :: error 00523 00524 CHARACTER(LEN=*), PARAMETER :: routineN = 'write_forces', 00525 routineP = moduleN//':'//routineN 00526 00527 CHARACTER(LEN=23) :: fmtstr3 00528 CHARACTER(LEN=36) :: fmtstr2 00529 CHARACTER(LEN=46) :: fmtstr1 00530 INTEGER :: i, ikind, iparticle, n 00531 LOGICAL :: failure 00532 00533 failure = .FALSE. 00534 00535 IF (output_unit > 0) THEN 00536 CPPrecondition(ASSOCIATED(particles),cp_failure_level,routineP,error,failure) 00537 n = MIN(MAX(1,ndigits),20) 00538 fmtstr1 = "(/,T2,A,/,/,T2,A,T11,A,T18,A,T35,A1,2( X,A1))" 00539 WRITE (UNIT=fmtstr1(39:40),FMT="(I2)") n + 6 00540 fmtstr2 = "(T2,I6,1X,I6,T21,A,T28,3(1X,F . ))" 00541 WRITE (UNIT=fmtstr2(33:34),FMT="(I2)") n 00542 WRITE (UNIT=fmtstr2(30:31),FMT="(I2)") n + 6 00543 fmtstr3 = "(T2,A,T28,4(1X,F . ))" 00544 WRITE (UNIT=fmtstr3(20:21),FMT="(I2)") n 00545 WRITE (UNIT=fmtstr3(17:18),FMT="(I2)") n + 6 00546 WRITE (UNIT=output_unit,FMT=fmtstr1)& 00547 label//" FORCES in [a.u.]","# Atom","Kind","Element","X","Y","Z" 00548 total_force(1:3) = 0.0_dp 00549 DO iparticle=1,particles%n_els 00550 ikind = particles%els(iparticle)%atomic_kind%kind_number 00551 IF (particles%els(iparticle)%atom_index /= 0) THEN 00552 i = particles%els(iparticle)%atom_index 00553 ELSE 00554 i = iparticle 00555 END IF 00556 WRITE (UNIT=output_unit,FMT=fmtstr2)& 00557 i,ikind,particles%els(iparticle)%atomic_kind%element_symbol,& 00558 particles%els(iparticle)%f(1:3) 00559 total_force(1:3) = total_force(1:3) + particles%els(iparticle)%f(1:3) 00560 END DO 00561 WRITE (UNIT=output_unit,FMT=fmtstr3)& 00562 "SUM OF "//label//" FORCES",total_force(1:3),SQRT(SUM(total_force(:)**2)) 00563 END IF 00564 00565 IF (PRESENT(grand_total_force)) THEN 00566 grand_total_force(1:3) = grand_total_force(1:3) + total_force(1:3) 00567 WRITE (UNIT=output_unit,FMT="(A)") "" 00568 WRITE (UNIT=output_unit,FMT=fmtstr3)& 00569 "GRAND TOTAL FORCE",grand_total_force(1:3),SQRT(SUM(grand_total_force(:)**2)) 00570 END IF 00571 00572 END SUBROUTINE write_forces 00573 00574 END MODULE force_env_utils
1.7.3