|
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 ! ***************************************************************************** 00013 MODULE reftraj_util 00014 00015 USE atomic_kind_list_types, ONLY: atomic_kind_list_type 00016 USE atomic_kind_types, ONLY: atomic_kind_type,& 00017 get_atomic_kind 00018 USE cp_files, ONLY: close_file,& 00019 open_file 00020 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 00021 cp_print_key_unit_nr 00022 USE cp_para_types, ONLY: cp_para_env_type 00023 USE cp_parser_methods, ONLY: parser_get_next_line 00024 USE cp_subsys_types, ONLY: cp_subsys_get,& 00025 cp_subsys_type 00026 USE cp_units, ONLY: cp_unit_to_cp2k 00027 USE distribution_1d_types, ONLY: distribution_1d_type 00028 USE f77_blas 00029 USE force_env_types, ONLY: force_env_get,& 00030 force_env_type 00031 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00032 section_vals_type,& 00033 section_vals_val_get 00034 USE kinds, ONLY: default_path_length,& 00035 default_string_length,& 00036 dp 00037 USE machine, ONLY: m_flush 00038 USE md_environment_types, ONLY: get_md_env,& 00039 md_environment_type 00040 USE message_passing, ONLY: mp_bcast,& 00041 mp_sum 00042 USE mol_kind_new_list_types, ONLY: mol_kind_new_list_type 00043 USE mol_new_list_types, ONLY: mol_new_list_type 00044 USE molecule_kind_types, ONLY: get_molecule_kind,& 00045 molecule_kind_type 00046 USE molecule_types_new, ONLY: get_molecule,& 00047 molecule_type 00048 USE particle_list_types, ONLY: particle_list_type 00049 USE particle_types, ONLY: particle_type 00050 USE physcon, ONLY: angstrom,& 00051 femtoseconds 00052 USE reftraj_types, ONLY: reftraj_msd_type,& 00053 reftraj_type 00054 USE simpar_types, ONLY: simpar_type 00055 USE termination, ONLY: stop_program 00056 USE util, ONLY: get_limit 00057 #include "cp_common_uses.h" 00058 00059 IMPLICIT NONE 00060 00061 PRIVATE 00062 00063 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'reftraj_util' 00064 00065 PUBLIC :: initialize_reftraj, compute_msd_reftraj, write_output_reftraj 00066 00067 CONTAINS 00068 00069 ! ***************************************************************************** 00076 SUBROUTINE initialize_reftraj(reftraj,reftraj_section,md_env,error) 00077 00078 TYPE(reftraj_type), POINTER :: reftraj 00079 TYPE(section_vals_type), POINTER :: reftraj_section 00080 TYPE(md_environment_type), POINTER :: md_env 00081 TYPE(cp_error_type), INTENT(inout) :: error 00082 00083 CHARACTER(len=*), PARAMETER :: routineN = 'initialize_reftraj', 00084 routineP = moduleN//':'//routineN 00085 00086 INTEGER :: natom, nline_to_skip, nskip 00087 LOGICAL :: failure, my_end 00088 TYPE(cp_para_env_type), POINTER :: para_env 00089 TYPE(cp_subsys_type), POINTER :: subsys 00090 TYPE(force_env_type), POINTER :: force_env 00091 TYPE(particle_list_type), POINTER :: particles 00092 TYPE(section_vals_type), POINTER :: msd_section 00093 TYPE(simpar_type), POINTER :: simpar 00094 00095 failure = .FALSE. 00096 00097 NULLIFY (force_env, msd_section, particles, simpar, subsys) 00098 CALL get_md_env(md_env=md_env, force_env=force_env, para_env=para_env,& 00099 simpar=simpar, error=error) 00100 CALL force_env_get(force_env=force_env, subsys=subsys,error=error) 00101 CALL cp_subsys_get(subsys=subsys, particles=particles,error=error) 00102 natom = particles%n_els 00103 00104 my_end = .FALSE. 00105 nline_to_skip = 0 00106 00107 nskip = reftraj%info%first_snapshot-1 00108 CPPostcondition(nskip>=0,cp_failure_level,routineP,error,failure) 00109 00110 IF(nskip > 0 ) THEN 00111 nline_to_skip = (natom+2)*nskip 00112 CALL parser_get_next_line(reftraj%info%traj_parser,nline_to_skip,at_end=my_end,error=error) 00113 END IF 00114 00115 reftraj%isnap = nskip 00116 CALL cp_assert(.NOT.my_end,cp_fatal_level,cp_assertion_failed,routineP,& 00117 "Reached the end of the trajectory file for REFTRAJ. Number of steps skipped "//& 00118 "equal to the number of steps present in the file. "//& 00119 CPSourceFileRef,& 00120 only_ionode=.TRUE.) 00121 00122 ! Cell File 00123 IF(reftraj%info%variable_volume) THEN 00124 IF(nskip > 0 ) THEN 00125 CALL parser_get_next_line(reftraj%info%cell_parser,nskip,at_end=my_end,error=error) 00126 END IF 00127 CALL cp_assert(.NOT.my_end,cp_fatal_level,cp_assertion_failed,routineP,& 00128 "Reached the end of the cell file for REFTRAJ. Number of steps skipped "//& 00129 "equal to the number of steps present in the file. "//& 00130 CPSourceFileRef,& 00131 only_ionode=.TRUE.) 00132 END IF 00133 00134 00135 reftraj%natom = natom 00136 IF(reftraj%info%last_snapshot>0) THEN 00137 simpar%nsteps = (reftraj%info%last_snapshot - reftraj%info%first_snapshot + 1) 00138 END IF 00139 00140 IF(reftraj%info%msd) THEN 00141 msd_section => section_vals_get_subs_vals(reftraj_section,"MSD",error=error) 00142 ! set up and printout 00143 CALL initialize_msd_reftraj(reftraj%msd,msd_section,reftraj,md_env,error=error) 00144 END IF 00145 00146 END SUBROUTINE initialize_reftraj 00147 00148 ! ***************************************************************************** 00155 SUBROUTINE initialize_msd_reftraj(msd,msd_section,reftraj,md_env,error) 00156 TYPE(reftraj_msd_type), POINTER :: msd 00157 TYPE(section_vals_type), POINTER :: msd_section 00158 TYPE(reftraj_type), POINTER :: reftraj 00159 TYPE(md_environment_type), POINTER :: md_env 00160 TYPE(cp_error_type), INTENT(inout) :: error 00161 00162 CHARACTER(len=*), PARAMETER :: routineN = 'initialize_msd_reftraj', 00163 routineP = moduleN//':'//routineN 00164 00165 CHARACTER(LEN=10) :: AA 00166 CHARACTER(LEN=80) :: title 00167 CHARACTER(LEN=default_path_length) :: filename 00168 CHARACTER(LEN=default_string_length) :: name 00169 INTEGER :: first_atom, iatom, ikind, imol, istat, last_atom, natom_read, 00170 nkind, nmol, nmolecule, nmolkind, npart 00171 LOGICAL :: failure 00172 REAL(KIND=dp) :: com(3), mass, mass_mol, tol, 00173 x, y, z 00174 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 00175 TYPE(cp_para_env_type), POINTER :: para_env 00176 TYPE(cp_subsys_type), POINTER :: subsys 00177 TYPE(force_env_type), POINTER :: force_env 00178 TYPE(mol_kind_new_list_type), POINTER :: molecule_kinds 00179 TYPE(mol_new_list_type), POINTER :: molecules 00180 TYPE(molecule_kind_type), DIMENSION(:), 00181 POINTER :: molecule_kind_set 00182 TYPE(molecule_kind_type), POINTER :: molecule_kind 00183 TYPE(molecule_type), DIMENSION(:), 00184 POINTER :: molecule_set 00185 TYPE(molecule_type), POINTER :: molecule 00186 TYPE(particle_list_type), POINTER :: particles 00187 TYPE(particle_type), DIMENSION(:), 00188 POINTER :: particle_set 00189 00190 failure = .FALSE. 00191 00192 NULLIFY (molecule, molecules, molecule_kind, molecule_kind_set,& 00193 molecule_kinds, molecule_set, subsys, force_env, particles, particle_set) 00194 CPPrecondition(.NOT. ASSOCIATED(msd),cp_failure_level,routineP,error,failure) 00195 00196 ALLOCATE(msd, STAT=istat) 00197 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00198 00199 NULLIFY(msd%ref0_pos) 00200 NULLIFY(msd%ref0_com_molecule) 00201 NULLIFY(msd%val_msd_kind) 00202 NULLIFY(msd%val_msd_molecule) 00203 NULLIFY(msd%disp_atom_index) 00204 NULLIFY(msd%disp_atom_dr) 00205 00206 CALL get_md_env(md_env=md_env, force_env=force_env, para_env=para_env,& 00207 error=error) 00208 CALL force_env_get(force_env=force_env,subsys=subsys,error=error) 00209 CALL cp_subsys_get(subsys=subsys, particles=particles,error=error) 00210 particle_set => particles%els 00211 npart = SIZE(particle_set,1) 00212 00213 msd%ref0_unit = -1 00214 CALL section_vals_val_get(msd_section,"REF0_FILENAME",c_val=filename,error=error) 00215 CALL open_file(TRIM(filename),unit_number=msd%ref0_unit) 00216 00217 ALLOCATE(msd%ref0_pos(3,reftraj%natom),STAT=istat) 00218 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00219 msd%ref0_pos = 0.0_dp 00220 00221 IF (para_env%mepos==para_env%source) THEN 00222 REWIND(msd%ref0_unit) 00223 READ(msd%ref0_unit,*,ERR=999,END=998) natom_read 00224 CPPostcondition(natom_read==reftraj%natom,cp_failure_level,routineP,error,failure) 00225 READ(msd%ref0_unit,'(A)',ERR=999,END=998) title 00226 msd%total_mass = 0.0_dp 00227 msd%ref0_com = 0.0_dp 00228 DO iatom = 1,natom_read 00229 READ(msd%ref0_unit,*,ERR=999,END=998) AA, x, y, z 00230 name=TRIM(particle_set(iatom)%atomic_kind%element_symbol) 00231 CPPostcondition((TRIM(AA)==name),cp_failure_level,routineP,error,failure) 00232 00233 x = cp_unit_to_cp2k(x,"angstrom",error=error) 00234 y = cp_unit_to_cp2k(y,"angstrom",error=error) 00235 z = cp_unit_to_cp2k(z,"angstrom",error=error) 00236 msd%ref0_pos(1,iatom) = x 00237 msd%ref0_pos(2,iatom) = y 00238 msd%ref0_pos(3,iatom) = z 00239 mass = particle_set(iatom)%atomic_kind%mass 00240 msd%ref0_com(1) = msd%ref0_com(1) + x * mass 00241 msd%ref0_com(2) = msd%ref0_com(2) + y * mass 00242 msd%ref0_com(3) = msd%ref0_com(3) + z * mass 00243 msd%total_mass = msd%total_mass + mass 00244 END DO 00245 msd%ref0_com = msd%ref0_com / msd%total_mass 00246 END IF 00247 CALL close_file(unit_number=msd%ref0_unit) 00248 00249 CALL mp_bcast(msd%total_mass,para_env%source,para_env%group) 00250 CALL mp_bcast(msd%ref0_pos,para_env%source,para_env%group) 00251 CALL mp_bcast(msd%ref0_com,para_env%source,para_env%group) 00252 00253 CALL section_vals_val_get(msd_section,"MSD_PER_KIND",l_val=msd%msd_kind,error=error) 00254 CALL section_vals_val_get(msd_section,"MSD_PER_MOLKIND",l_val=msd%msd_molecule,error=error) 00255 CALL section_vals_val_get(msd_section,"MSD_PER_REGION",l_val=msd%msd_region,error=error) 00256 00257 CALL section_vals_val_get(msd_section,"DISPLACED_ATOM",l_val=msd%disp_atom,error=error) 00258 IF(msd%disp_atom) THEN 00259 ALLOCATE(msd%disp_atom_index(npart),STAT=istat) 00260 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00261 msd%disp_atom_index = 0 00262 ALLOCATE(msd%disp_atom_dr(3,npart),STAT=istat) 00263 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00264 msd%disp_atom_dr = 0.0_dp 00265 msd%msd_kind = .TRUE. 00266 END IF 00267 CALL section_vals_val_get(msd_section,"DISPLACEMENT_TOL",r_val=tol,error=error) 00268 msd%disp_atom_tol = tol*tol 00269 00270 IF(msd%msd_kind) THEN 00271 CALL cp_subsys_get(subsys=subsys,atomic_kinds=atomic_kinds, error=error) 00272 nkind = atomic_kinds%n_els 00273 00274 ALLOCATE(msd%val_msd_kind(4,nkind),STAT=istat) 00275 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00276 msd%val_msd_kind = 0.0_dp 00277 END IF 00278 00279 IF(msd%msd_molecule) THEN 00280 CALL cp_subsys_get(subsys=subsys, molecules_new=molecules,& 00281 molecule_kinds_new=molecule_kinds,error=error) 00282 nmolkind = molecule_kinds%n_els 00283 ALLOCATE(msd%val_msd_molecule(4,nmolkind),STAT=istat) 00284 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00285 00286 molecule_kind_set => molecule_kinds%els 00287 molecule_set => molecules%els 00288 nmol = molecules%n_els 00289 00290 ALLOCATE(msd%ref0_com_molecule(3,nmol), STAT=istat) 00291 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00292 00293 DO ikind = 1, nmolkind 00294 molecule_kind => molecule_kind_set(ikind) 00295 CALL get_molecule_kind (molecule_kind=molecule_kind, nmolecule=nmolecule ) 00296 DO imol = 1,nmolecule 00297 molecule => molecule_set(molecule_kind%molecule_list(imol)) 00298 CALL get_molecule ( molecule=molecule, first_atom = first_atom, last_atom = last_atom ) 00299 com = 0.0_dp 00300 mass_mol = 0.0_dp 00301 DO iatom = first_atom, last_atom 00302 mass = particle_set(iatom)%atomic_kind%mass 00303 com(1) = com(1) + msd%ref0_pos(1,iatom)*mass 00304 com(2) = com(2) + msd%ref0_pos(2,iatom)*mass 00305 com(3) = com(3) + msd%ref0_pos(3,iatom)*mass 00306 mass_mol = mass_mol+mass 00307 ENDDO ! iatom 00308 msd%ref0_com_molecule(1,molecule_kind%molecule_list(imol)) = com(1)/mass_mol 00309 msd%ref0_com_molecule(2,molecule_kind%molecule_list(imol)) = com(2)/mass_mol 00310 msd%ref0_com_molecule(3,molecule_kind%molecule_list(imol)) = com(3)/mass_mol 00311 END DO ! imol 00312 ENDDO ! ikind 00313 END IF 00314 00315 IF(msd%msd_region) THEN 00316 00317 END IF 00318 00319 RETURN 00320 998 CONTINUE ! end of file 00321 CALL stop_program(routineN,moduleN,__LINE__,& 00322 "End of reference positions file reached") 00323 999 CONTINUE ! error 00324 CALL stop_program(routineN,moduleN,__LINE__,& 00325 "Error reading reference positions file") 00326 00327 END SUBROUTINE initialize_msd_reftraj 00328 00329 ! ***************************************************************************** 00336 SUBROUTINE compute_msd_reftraj(reftraj,md_env,particle_set,error) 00337 00338 TYPE(reftraj_type), POINTER :: reftraj 00339 TYPE(md_environment_type), POINTER :: md_env 00340 TYPE(particle_type), DIMENSION(:), 00341 POINTER :: particle_set 00342 TYPE(cp_error_type), INTENT(inout) :: error 00343 00344 CHARACTER(len=*), PARAMETER :: routineN = 'compute_msd_reftraj', 00345 routineP = moduleN//':'//routineN 00346 00347 INTEGER :: atom, bo(2), first_atom, iatom, ikind, imol, imol_global, 00348 last_atom, mepos, natom_kind, nmol_per_kind, nmolecule, nmolkind, num_pe 00349 INTEGER, DIMENSION(:), POINTER :: atom_list 00350 LOGICAL :: failure 00351 REAL(KIND=dp) :: com(3), diff2_com(4), dr2, 00352 dx, dy, dz, mass, mass_mol, 00353 msd_mkind(4), rcom(3) 00354 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 00355 TYPE(atomic_kind_type), POINTER :: atomic_kind 00356 TYPE(cp_para_env_type), POINTER :: para_env 00357 TYPE(cp_subsys_type), POINTER :: subsys 00358 TYPE(distribution_1d_type), POINTER :: local_molecules 00359 TYPE(force_env_type), POINTER :: force_env 00360 TYPE(mol_kind_new_list_type), POINTER :: molecule_kinds 00361 TYPE(mol_new_list_type), POINTER :: molecules 00362 TYPE(molecule_kind_type), DIMENSION(:), 00363 POINTER :: molecule_kind_set 00364 TYPE(molecule_kind_type), POINTER :: molecule_kind 00365 TYPE(molecule_type), DIMENSION(:), 00366 POINTER :: molecule_set 00367 TYPE(molecule_type), POINTER :: molecule 00368 00369 failure = .FALSE. 00370 NULLIFY(force_env,para_env,subsys) 00371 NULLIFY(atomic_kind,atomic_kinds,atom_list) 00372 NULLIFY(local_molecules, molecule, molecule_kind, molecule_kinds,& 00373 molecule_kind_set, molecules, molecule_set) 00374 00375 CALL get_md_env(md_env=md_env, force_env=force_env, para_env=para_env,& 00376 error=error) 00377 CALL force_env_get(force_env=force_env,subsys=subsys,& 00378 error=error) 00379 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds,error=error) 00380 00381 num_pe = para_env%num_pe 00382 mepos = para_env%mepos 00383 00384 00385 IF(reftraj%msd%msd_kind)THEN 00386 reftraj%msd%val_msd_kind = 0.0_dp 00387 reftraj%msd%num_disp_atom = 0 00388 reftraj%msd%disp_atom_dr = 0.0_dp 00389 ! compute com 00390 rcom = 0.0_dp 00391 DO ikind=1,atomic_kinds%n_els 00392 atomic_kind => atomic_kinds%els(ikind) 00393 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00394 atom_list=atom_list,& 00395 natom=natom_kind, mass=mass) 00396 bo = get_limit(natom_kind,num_pe,mepos) 00397 DO iatom = bo(1), bo(2) 00398 atom = atom_list(iatom) 00399 rcom(1) = rcom(1) + particle_set(atom)%r(1)*mass 00400 rcom(2) = rcom(2) + particle_set(atom)%r(2)*mass 00401 rcom(3) = rcom(3) + particle_set(atom)%r(3)*mass 00402 END DO 00403 END DO 00404 CALL mp_sum(rcom,para_env%group) 00405 rcom = rcom / reftraj%msd%total_mass 00406 reftraj%msd%drcom (1) = rcom(1) - reftraj%msd%ref0_com(1) 00407 reftraj%msd%drcom (2) = rcom(2) - reftraj%msd%ref0_com(2) 00408 reftraj%msd%drcom (3) = rcom(3) - reftraj%msd%ref0_com(3) 00409 ! IF(para_env%ionode) WRITE(6,'(A,T50,3f10.5)') ' COM displacement (dx,dy,dz) [angstrom]: ', & 00410 ! drcom(1)*angstrom,drcom(2)*angstrom,drcom(3)*angstrom 00411 ! compute_com 00412 00413 DO ikind=1,atomic_kinds%n_els 00414 atomic_kind => atomic_kinds%els(ikind) 00415 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00416 atom_list=atom_list,& 00417 natom=natom_kind) 00418 bo = get_limit(natom_kind,num_pe,mepos) 00419 DO iatom = bo(1), bo(2) 00420 atom = atom_list(iatom) 00421 dx = particle_set(atom)%r(1)-reftraj%msd%ref0_pos(1,atom) - & 00422 reftraj%msd%drcom(1) 00423 dy = particle_set(atom)%r(2)-reftraj%msd%ref0_pos(2,atom) - & 00424 reftraj%msd%drcom(2) 00425 dz = particle_set(atom)%r(3)-reftraj%msd%ref0_pos(3,atom) - & 00426 reftraj%msd%drcom(3) 00427 dr2 = dx*dx + dy*dy + dz*dz 00428 00429 reftraj%msd%val_msd_kind(1,ikind) = reftraj%msd%val_msd_kind(1,ikind) + dx*dx 00430 reftraj%msd%val_msd_kind(2,ikind) = reftraj%msd%val_msd_kind(2,ikind) + dy*dy 00431 reftraj%msd%val_msd_kind(3,ikind) = reftraj%msd%val_msd_kind(3,ikind) + dz*dz 00432 reftraj%msd%val_msd_kind(4,ikind) = reftraj%msd%val_msd_kind(4,ikind) + dr2 00433 00434 IF(reftraj%msd%disp_atom) THEN 00435 IF( dr2 > reftraj%msd%disp_atom_tol ) THEN 00436 reftraj%msd%num_disp_atom = reftraj%msd%num_disp_atom +1 00437 reftraj%msd%disp_atom_dr(1, atom) = dx 00438 reftraj%msd%disp_atom_dr(2, atom) = dy 00439 reftraj%msd%disp_atom_dr(3, atom) = dz 00440 END IF 00441 END IF 00442 END DO !iatom 00443 reftraj%msd%val_msd_kind(1:4,ikind) = & 00444 reftraj%msd%val_msd_kind(1:4,ikind)/REAL(natom_kind,KIND=dp) 00445 00446 END DO ! ikind 00447 ENDIF 00448 CALL mp_sum(reftraj%msd%val_msd_kind,para_env%group) 00449 CALL mp_sum(reftraj%msd%num_disp_atom,para_env%group) 00450 CALL mp_sum(reftraj%msd%disp_atom_dr,para_env%group) 00451 00452 IF(reftraj%msd%msd_molecule) THEN 00453 CALL cp_subsys_get(subsys=subsys, local_molecules_new=local_molecules, & 00454 molecules_new=molecules, molecule_kinds_new=molecule_kinds, error=error) 00455 00456 nmolkind = molecule_kinds%n_els 00457 molecule_kind_set => molecule_kinds%els 00458 molecule_set => molecules%els 00459 00460 reftraj%msd%val_msd_molecule = 0.0_dp 00461 DO ikind = 1,nmolkind 00462 molecule_kind => molecule_kind_set(ikind) 00463 CALL get_molecule_kind (molecule_kind=molecule_kind, nmolecule=nmolecule ) 00464 nmol_per_kind = local_molecules%n_el(ikind) 00465 msd_mkind = 0.0_dp 00466 DO imol = 1, nmol_per_kind 00467 imol_global = local_molecules%list(ikind)%array(imol) 00468 molecule => molecule_set ( imol_global ) 00469 CALL get_molecule (molecule,first_atom=first_atom,last_atom=last_atom) 00470 00471 com = 0.0_dp 00472 mass_mol = 0.0_dp 00473 DO iatom = first_atom, last_atom 00474 mass = particle_set(iatom)%atomic_kind%mass 00475 com(1) = com(1) + particle_set(iatom)%r(1)*mass 00476 com(2) = com(2) + particle_set(iatom)%r(2)*mass 00477 com(3) = com(3) + particle_set(iatom)%r(3)*mass 00478 mass_mol = mass_mol+mass 00479 ENDDO ! iatom 00480 com(1) = com(1)/mass_mol 00481 com(2) = com(2)/mass_mol 00482 com(3) = com(3)/mass_mol 00483 diff2_com(1) = com(1)- reftraj%msd%ref0_com_molecule(1,imol_global) 00484 diff2_com(2) = com(2)- reftraj%msd%ref0_com_molecule(2,imol_global) 00485 diff2_com(3) = com(3)- reftraj%msd%ref0_com_molecule(3,imol_global) 00486 diff2_com(1) = diff2_com(1)*diff2_com(1) 00487 diff2_com(2) = diff2_com(2)*diff2_com(2) 00488 diff2_com(3) = diff2_com(3)*diff2_com(3) 00489 diff2_com(4) = diff2_com(1) + diff2_com(2) + diff2_com(3) 00490 msd_mkind(1) = msd_mkind(1) + diff2_com(1) 00491 msd_mkind(2) = msd_mkind(2) + diff2_com(2) 00492 msd_mkind(3) = msd_mkind(3) + diff2_com(3) 00493 msd_mkind(4) = msd_mkind(4) + diff2_com(4) 00494 ENDDO ! imol 00495 00496 reftraj%msd%val_msd_molecule(1,ikind) = msd_mkind(1)/REAL(nmolecule,KIND=dp) 00497 reftraj%msd%val_msd_molecule(2,ikind) = msd_mkind(2)/REAL(nmolecule,KIND=dp) 00498 reftraj%msd%val_msd_molecule(3,ikind) = msd_mkind(3)/REAL(nmolecule,KIND=dp) 00499 reftraj%msd%val_msd_molecule(4,ikind) = msd_mkind(4)/REAL(nmolecule,KIND=dp) 00500 END DO ! ikind 00501 CALL mp_sum(reftraj%msd%val_msd_molecule, para_env%group) 00502 00503 END IF 00504 00505 END SUBROUTINE compute_msd_reftraj 00506 00507 ! ***************************************************************************** 00514 SUBROUTINE write_output_reftraj(md_env,error) 00515 TYPE(md_environment_type), POINTER :: md_env 00516 TYPE(cp_error_type), INTENT(inout) :: error 00517 00518 CHARACTER(len=*), PARAMETER :: routineN = 'write_output_reftraj', 00519 routineP = moduleN//':'//routineN 00520 00521 CHARACTER(LEN=default_string_length) :: my_act, my_mittle, my_pos 00522 INTEGER :: iat, ikind, nkind, out_msd 00523 LOGICAL, SAVE :: first_entry = .FALSE. 00524 TYPE(cp_logger_type), POINTER :: logger 00525 TYPE(force_env_type), POINTER :: force_env 00526 TYPE(reftraj_type), POINTER :: reftraj 00527 TYPE(section_vals_type), POINTER :: reftraj_section, root_section 00528 00529 NULLIFY(logger) 00530 logger => cp_error_get_logger(error) 00531 00532 NULLIFY(reftraj) 00533 NULLIFY(reftraj_section,root_section) 00534 00535 CALL get_md_env(md_env=md_env, force_env=force_env, & 00536 reftraj=reftraj, error=error) 00537 00538 CALL force_env_get(force_env=force_env,root_section=root_section,& 00539 error=error) 00540 00541 reftraj_section => section_vals_get_subs_vals(root_section,& 00542 "MOTION%MD%REFTRAJ",error=error) 00543 00544 my_pos = "APPEND" 00545 my_act = "WRITE" 00546 00547 IF(reftraj%init.AND.(reftraj%isnap==reftraj%info%first_snapshot)) THEN 00548 my_pos = "REWIND" 00549 first_entry = .TRUE. 00550 END IF 00551 00552 IF(reftraj%info%msd) THEN 00553 IF(reftraj%msd%msd_kind)THEN 00554 nkind = SIZE(reftraj%msd%val_msd_kind,2) 00555 DO ikind = 1,nkind 00556 my_mittle ="k"//TRIM(ADJUSTL(cp_to_string(ikind))) 00557 out_msd = cp_print_key_unit_nr(logger,reftraj_section,"PRINT%MSD_KIND",& 00558 extension=".msd", file_position=my_pos, file_action=my_act,& 00559 file_form="FORMATTED", middle_name=TRIM(my_mittle), & 00560 error=error) 00561 IF(out_msd>0) THEN 00562 WRITE(UNIT=out_msd,FMT="(I8, F12.3,4F20.10)") reftraj%itimes, & 00563 reftraj%time*femtoseconds, & 00564 reftraj%msd%val_msd_kind(1:4,ikind)*angstrom*angstrom 00565 CALL m_flush(out_msd) 00566 END IF 00567 CALL cp_print_key_finished_output(out_msd,logger,reftraj_section,& 00568 "PRINT%MSD_KIND", error=error) 00569 END DO 00570 END IF 00571 IF(reftraj%msd%msd_molecule) THEN 00572 nkind = SIZE(reftraj%msd%val_msd_molecule,2) 00573 DO ikind = 1,nkind 00574 my_mittle ="mk"//TRIM(ADJUSTL(cp_to_string(ikind))) 00575 out_msd = cp_print_key_unit_nr(logger,reftraj_section,"PRINT%MSD_MOLECULE",& 00576 extension=".msd", file_position=my_pos, file_action=my_act,& 00577 file_form="FORMATTED", middle_name=TRIM(my_mittle), & 00578 error=error) 00579 IF(out_msd>0) THEN 00580 WRITE(UNIT=out_msd,FMT="(I8, F12.3,4F20.10)") reftraj%itimes, & 00581 reftraj%time*femtoseconds, & 00582 reftraj%msd%val_msd_molecule(1:4,ikind)*angstrom*angstrom 00583 CALL m_flush(out_msd) 00584 END IF 00585 CALL cp_print_key_finished_output(out_msd,logger,reftraj_section,& 00586 "PRINT%MSD_MOLECULE", error=error) 00587 END DO 00588 END IF 00589 IF(reftraj%msd%disp_atom) THEN 00590 00591 IF(first_entry) my_pos = "REWIND" 00592 my_mittle ="disp_at" 00593 out_msd = cp_print_key_unit_nr(logger,reftraj_section,"PRINT%DISPLACED_ATOM",& 00594 extension=".msd", file_position=my_pos, file_action=my_act,& 00595 file_form="FORMATTED", middle_name=TRIM(my_mittle), & 00596 error=error) 00597 IF(out_msd>0 .AND. reftraj%msd%num_disp_atom>0) THEN 00598 IF(first_entry) THEN 00599 first_entry = .FALSE. 00600 END IF 00601 WRITE(UNIT=out_msd,FMT="(A,T7,I8, A, T29, F12.3, A, T50, I10)") "# i = ", reftraj%itimes, " time (fs) = ",& 00602 reftraj%time*femtoseconds, " nat = ", reftraj%msd%num_disp_atom 00603 DO iat = 1,SIZE(reftraj%msd%disp_atom_dr,2) 00604 IF(ABS(reftraj%msd%disp_atom_dr(1,iat)) > 0.0_dp) THEN 00605 WRITE(UNIT=out_msd,FMT="(I8, 3F20.10)") iat,& !reftraj%msd%disp_atom_index(iat),& 00606 reftraj%msd%disp_atom_dr(1,iat)*angstrom, & 00607 reftraj%msd%disp_atom_dr(2,iat)*angstrom,& 00608 reftraj%msd%disp_atom_dr(3,iat)*angstrom 00609 END IF 00610 END DO 00611 ENDIF 00612 CALL cp_print_key_finished_output(out_msd,logger,reftraj_section,& 00613 "PRINT%DISPLACED_ATOM", error=error) 00614 END IF 00615 ENDIF ! msd 00616 reftraj%init = .FALSE. 00617 00618 END SUBROUTINE write_output_reftraj 00619 00620 END MODULE reftraj_util 00621
1.7.3