diff --git a/src/fist_nonbond_force.F b/src/fist_nonbond_force.F index 067bffae5a..92bfeb5f9c 100644 --- a/src/fist_nonbond_force.F +++ b/src/fist_nonbond_force.F @@ -96,7 +96,7 @@ SUBROUTINE force_nonbond ( ewald_param, part, pnode, box, potparm, & !MK TYPE(neighbor_list_set_type), POINTER :: neighbor_list_set - INTEGER :: ilist,inode,nlist,nnode + INTEGER :: ilist,nlist,nnode REAL(dbl), DIMENSION(3) :: dra,drb,rab_last_update diff --git a/src/force_control.F b/src/force_control.F index 44f1c702ec..5a410c7189 100644 --- a/src/force_control.F +++ b/src/force_control.F @@ -30,6 +30,7 @@ MODULE force_control USE kinds, ONLY : dbl USE qs_build_ks_matrix, ONLY : qs_ks_env_type USE qs_environment_types, ONLY : qs_environment_type + USE qs_force, ONLY: qs_forces_update USE termination, ONLY : stop_program IMPLICIT NONE @@ -79,6 +80,7 @@ SUBROUTINE force ( globenv, box_change ) INTEGER :: beads, i !------------------------------------------------------------------------------ + SELECT CASE ( in_use ) CASE ( FIST ) SELECT CASE ( fist_env % pimd ) @@ -100,6 +102,7 @@ SUBROUTINE force ( globenv, box_change ) box_change, globenv ) END SELECT CASE ( QUICKSTEP ) + CALL qs_forces_update(qs_env,globenv,ks_env) CASE ( POL ) CASE ( KIM_GORDON ) CALL kg_force_control ( kg_env, globenv ) diff --git a/src/particle_types.F b/src/particle_types.F index d1358ca226..e8c8f190d8 100644 --- a/src/particle_types.F +++ b/src/particle_types.F @@ -221,7 +221,9 @@ CONTAINS ! *** Local variables *** TYPE(gto_basis_set_type), POINTER :: orb_basis_set - INTEGER :: iatom,icgf,isgf,istat,natom,ncgf,nsgf + + REAL(wp) :: mass + INTEGER :: iatom,icgf,isgf,istat,kind_number,natom,ncgf,nsgf ! --------------------------------------------------------------------------- @@ -231,7 +233,14 @@ CONTAINS isgf = 0 DO iatom=1,natom CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& - orb_basis_set=orb_basis_set) + orb_basis_set=orb_basis_set,& + kind_number=kind_number,& + mass=mass) + ALLOCATE (particle_set(iatom)%prop,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"particle_set(iatom)%prop",0) + particle_set(iatom)%prop%charge = 0.0_wp + particle_set(iatom)%prop%mass = mass + particle_set(iatom)%prop%ptype = kind_number CALL get_gto_basis_set(gto_basis_set=orb_basis_set,ncgf=ncgf,nsgf=nsgf) particle_set(iatom)%first_cgf = icgf + 1 icgf = icgf + ncgf @@ -570,7 +579,7 @@ CONTAINS rab_image(:) = rb_image(:) - ra(:) dab_image = SQRT(rab_image(1)*rab_image(1) +& rab_image(2)*rab_image(2) +& - rab_image(3)*rab_image(3)) + rab_image(3)*rab_image(3)) IF (dab_image < dab) dab = dab_image END DO END DO diff --git a/src/qs_energy.F b/src/qs_energy.F index 5c738c8315..a5aa0a2511 100644 --- a/src/qs_energy.F +++ b/src/qs_energy.F @@ -49,7 +49,8 @@ MODULE qs_energy CHARACTER(LEN=*), PARAMETER :: module_name = "qs_energy" - PUBLIC :: qs_energies + PUBLIC :: qs_energies,& + qs_energies_update ! ***************************************************************************** @@ -65,9 +66,9 @@ CONTAINS ! *************************************************************************** - TYPE(global_environment_type), INTENT(INOUT) :: globenv - TYPE(qs_environment_type), INTENT(INOUT) :: qs_env - TYPE(qs_ks_env_type), INTENT(INOUT) :: ks_env + TYPE(global_environment_type), INTENT(IN) :: globenv + TYPE(qs_environment_type), INTENT(INOUT) :: qs_env + TYPE(qs_ks_env_type), INTENT(INOUT) :: ks_env ! *** Local parameters *** @@ -130,9 +131,9 @@ CONTAINS ! *************************************************************************** - TYPE(global_environment_type), INTENT(INOUT) :: globenv - TYPE(qs_environment_type), INTENT(INOUT) :: qs_env - TYPE(qs_ks_env_type), INTENT(INOUT) :: ks_env + TYPE(global_environment_type), INTENT(IN) :: globenv + TYPE(qs_environment_type), INTENT(INOUT) :: qs_env + TYPE(qs_ks_env_type), INTENT(INOUT) :: ks_env ! *** Local parameters *** @@ -167,10 +168,6 @@ CONTAINS globenv=globenv,& calculate_forces=.FALSE.) -! *** Updates rho structure *** - - CALL qs_env_rebuild_rho(qs_env=qs_env) - ! *** Perform a SCF run *** CALL scf(ks_env,qs_env,globenv) diff --git a/src/qs_force.F b/src/qs_force.F index ed64bda13d..444f22bbdc 100644 --- a/src/qs_force.F +++ b/src/qs_force.F @@ -28,12 +28,14 @@ MODULE qs_force wp => dp,& wp_size => dp_size USE message_passing, ONLY: mp_sum + USE particle_types, ONLY: particle_type USE qs_build_ks_matrix, ONLY: qs_ks_did_change,& qs_ks_env_type,& qs_ks_update_qs_env USE qs_core_energies, ONLY: calculate_ecore_overlap USE qs_core_hamiltonian, ONLY: build_core_hamiltonian_matrix - USE qs_energy, ONLY: qs_energies + USE qs_energy, ONLY: qs_energies,& + qs_energies_update USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type USE qs_force_types, ONLY: init_qs_force,& @@ -64,7 +66,8 @@ MODULE qs_force ! *** Public subroutines *** - PUBLIC :: qs_forces + PUBLIC :: qs_forces,& + qs_forces_update ! ***************************************************************************** @@ -80,9 +83,9 @@ CONTAINS ! *************************************************************************** - TYPE(global_environment_type), INTENT(INOUT) :: globenv - TYPE(qs_environment_type), INTENT(INOUT) :: qs_env - TYPE(qs_ks_env_type), INTENT(INOUT) :: ks_env + TYPE(global_environment_type), INTENT(IN) :: globenv + TYPE(qs_environment_type), INTENT(INOUT) :: qs_env + TYPE(qs_ks_env_type), INTENT(INOUT) :: ks_env ! *** Local parameters *** @@ -92,16 +95,19 @@ CONTAINS ! *** Local variables *** - TYPE(mo_set_type), POINTER :: mo_set - TYPE(real_matrix_type), POINTER :: w + TYPE(mo_set_type), POINTER :: mo_set + TYPE(qs_rho_type), POINTER :: rho + TYPE(real_matrix_type), POINTER :: w - INTEGER :: handle,ikind,istat,output_unit + INTEGER :: handle,i,iatom,ikind,istat,natom,output_unit TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(mo_set_p_type), DIMENSION(:), POINTER :: c + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(qs_force_type), DIMENSION(:), POINTER :: force TYPE(real_matrix_p_type), DIMENSION(:), POINTER :: k,p,s - type(qs_rho_type), pointer :: rho + + INTEGER, DIMENSION(:), ALLOCATABLE :: atom_of_kind,kind_of ! --------------------------------------------------------------------------- @@ -113,12 +119,23 @@ CONTAINS CALL get_qs_env(qs_env=qs_env,& atomic_kind_set=atomic_kind_set,& + particle_set=particle_set,& s=s,& c=c,& k=k,& rho=rho,& force=force) + natom = SIZE(particle_set) + + ALLOCATE (atom_of_kind(natom),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"atom_of_kind",natom*int_size) + ALLOCATE (kind_of(natom),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"kind_of",natom*int_size) + CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& + atom_of_kind=atom_of_kind,& + kind_of=kind_of) + p => rho%rho_ao CALL init_qs_force(force) @@ -169,16 +186,163 @@ CONTAINS force(ikind)%rho_elec(:,:) END DO + DO iatom=1,natom + ikind = kind_of(iatom) + i = atom_of_kind(iatom) + particle_set(iatom)%f(1:3) = force(ikind)%total(1:3,i) + END DO + CALL write_forces(force,atomic_kind_set,globenv) CALL deallocate_matrix(w) + DEALLOCATE (atom_of_kind,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"atom_of_kind") + + DEALLOCATE (kind_of,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"kind_of") + CALL timestop(0.0_wp,handle) CALL write_checkpoint_information("leaving "//routine_name,globenv) END SUBROUTINE qs_forces +! ***************************************************************************** + + SUBROUTINE qs_forces_update(qs_env,globenv,ks_env) + +! Purpose: Calculate the Quickstep forces. + +! History: - Creation (29.10.2002,MK) + +! *************************************************************************** + + TYPE(global_environment_type), INTENT(IN) :: globenv + TYPE(qs_environment_type), INTENT(INOUT) :: qs_env + TYPE(qs_ks_env_type), INTENT(INOUT) :: ks_env + +! *** Local parameters *** + + CHARACTER(LEN=*), PARAMETER :: routine_name = "qs_forces" + CHARACTER(LEN=*), PARAMETER :: routine =& + "SUBROUTINE "//routine_name//" (MODULE "//module_name//")" + +! *** Local variables *** + + TYPE(mo_set_type), POINTER :: mo_set + TYPE(qs_rho_type), POINTER :: rho + TYPE(real_matrix_type), POINTER :: w + + INTEGER :: handle,i,iatom,ikind,istat,natom,output_unit + + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(mo_set_p_type), DIMENSION(:), POINTER :: c + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(qs_force_type), DIMENSION(:), POINTER :: force + TYPE(real_matrix_p_type), DIMENSION(:), POINTER :: k,p,s + + INTEGER, DIMENSION(:), ALLOCATABLE :: atom_of_kind,kind_of + +! --------------------------------------------------------------------------- + + CALL write_checkpoint_information("entering "//routine_name,globenv) + + CALL timeset(routine_name,"I","",handle) + + CALL qs_energies_update(qs_env,globenv,ks_env) + + CALL get_qs_env(qs_env=qs_env,& + atomic_kind_set=atomic_kind_set,& + particle_set=particle_set,& + s=s,& + c=c,& + k=k,& + rho=rho,& + force=force) + + natom = SIZE(particle_set) + + ALLOCATE (atom_of_kind(natom),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"atom_of_kind",natom*int_size) + ALLOCATE (kind_of(natom),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"kind_of",natom*int_size) + CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& + atom_of_kind=atom_of_kind,& + kind_of=kind_of) + + p => rho%rho_ao + + CALL init_qs_force(force) + +! *** Force component from the core charge distribution overlap *** + + CALL calculate_ecore_overlap(qs_env=qs_env,& + globenv=globenv,& + calculate_forces=.TRUE.) + +! *** Load MO set (QS wavefunction) *** + + mo_set => c(1)%mo_set + + NULLIFY (w) + CALL replicate_matrix_structure(s(1)%matrix,w,"W MATRIX") + + CALL calculate_w_matrix(mo_set,k(1)%matrix,w,globenv) + IF (globenv%ionode.AND.globenv%print%w_matrix) THEN + CALL write_sparse_matrix(w,4,6,qs_env,globenv) + END IF + + CALL build_core_hamiltonian_matrix(qs_env=qs_env,& + globenv=globenv,& + calculate_forces=.TRUE.,& + p=p(1)%matrix,& + w=w) + CALL qs_ks_did_change(ks_env,qs_env,globenv,p_changed=.true.) + CALL qs_ks_update_qs_env(ks_env=ks_env,& + qs_env=qs_env,& + global_env=globenv,& + rebuild=.TRUE.,& + calculate_forces=.TRUE.) + + DO ikind=1,SIZE(force) + CALL mp_sum(force(ikind)%gth_ppl,globenv%group) + CALL mp_sum(force(ikind)%gth_ppnl,globenv%group) + CALL mp_sum(force(ikind)%kinetic,globenv%group) + CALL mp_sum(force(ikind)%overlap,globenv%group) + CALL mp_sum(force(ikind)%rho_elec,globenv%group) + force(ikind)%total(:,:) = force(ikind)%total(:,:) +& + force(ikind)%core_overlap(:,:) +& + force(ikind)%gth_ppl(:,:) +& + force(ikind)%gth_ppnl(:,:) +& + force(ikind)%kinetic(:,:) +& + force(ikind)%overlap(:,:) +& + force(ikind)%rho_core(:,:) +& + force(ikind)%rho_elec(:,:) + END DO + + DO iatom=1,natom + ikind = kind_of(iatom) + i = atom_of_kind(iatom) + particle_set(iatom)%f(1:3) = force(ikind)%total(1:3,i) + END DO + + CALL write_forces(force,atomic_kind_set,globenv) + + CALL deallocate_matrix(w) + + DEALLOCATE (atom_of_kind,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"atom_of_kind") + + DEALLOCATE (kind_of,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"kind_of") + + CALL timestop(0.0_wp,handle) + + CALL write_checkpoint_information("leaving "//routine_name,globenv) + + END SUBROUTINE qs_forces_update + ! ***************************************************************************** SUBROUTINE write_forces(qs_force,atomic_kind_set,globenv) @@ -202,8 +366,9 @@ CONTAINS INTEGER :: i,iatom,ikind,istat,natom,output_unit + REAL(wp), DIMENSION(3) :: grand_total + INTEGER, DIMENSION(:), ALLOCATABLE :: atom_of_kind,kind_of - REAL(wp) :: grand_total(3) ! --------------------------------------------------------------------------- @@ -216,7 +381,6 @@ CONTAINS "and cannot be printed") END IF - output_unit = globenv%scr output_unit = globenv%scr CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& @@ -233,7 +397,8 @@ CONTAINS FMT="(/,/,T2,A,/,/,T3,A,T11,A,T23,A,T43,A1,2(13X,A1))")& "FORCES","Atom","Kind","Component","X","Y","Z" - grand_total(:)=0.0_wp + grand_total(:) = 0.0_wp + DO iatom=1,natom ikind = kind_of(iatom) i = atom_of_kind(iatom) @@ -246,10 +411,17 @@ CONTAINS iatom,ikind," rho_core",qs_force(ikind)%rho_core(1:3,i),& iatom,ikind," rho_elec",qs_force(ikind)%rho_elec(1:3,i),& iatom,ikind," total",qs_force(ikind)%total(1:3,i) - grand_total(1:3)=grand_total(1:3)+qs_force(ikind)%total(1:3,i) + grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3,i) END DO - WRITE(UNIT=output_unit,FMT="(/,(T2,A,T34,3F14.6))") "Sum of total", & - grand_total(1:3) + + WRITE (UNIT=output_unit,FMT="(/,(T2,A,T34,3F14.6))")& + "Sum of total",grand_total(1:3) + + DEALLOCATE (atom_of_kind,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"atom_of_kind") + + DEALLOCATE (kind_of,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"kind_of") END SUBROUTINE write_forces diff --git a/src/qs_main.F b/src/qs_main.F index bdb4e83cd6..53543852ec 100644 --- a/src/qs_main.F +++ b/src/qs_main.F @@ -29,7 +29,7 @@ MODULE qs_main USE checkpoint_handler, ONLY: write_checkpoint_information USE global_types, ONLY: global_environment_type - USE kinds, ONLY: wp => dp + USE kinds, ONLY: dp USE qs_md, ONLY: qs_mol_dyn USE qs_build_ks_matrix, ONLY: qs_ks_env_type USE qs_energy, ONLY: qs_energies @@ -110,7 +110,7 @@ CONTAINS CALL qs_finalize(qs_env,globenv) - CALL timestop(0.0_wp,handle) + CALL timestop(0.0_dp,handle) CALL write_checkpoint_information("leaving "//routine_name,globenv) diff --git a/src/qs_md.F b/src/qs_md.F index 57d3a6a711..169facfd9b 100644 --- a/src/qs_md.F +++ b/src/qs_md.F @@ -21,16 +21,58 @@ MODULE qs_md + USE atoms_input, ONLY: set_particle_properties USE checkpoint_handler, ONLY: write_checkpoint_information + USE cntl_input, ONLY: read_cntl_section + USE convert_units, ONLY: convert + USE extended_system_types, ONLY: extended_system_type + USE force_control, ONLY: force,& + set_force_env USE global_types, ONLY: global_environment_type - USE kinds, ONLY: wp => dp + USE initialize_extended_types, ONLY: initialize_nhc_part,& + initialize_nhc_baro,& + initialize_npt_type + USE initialize_molecule_types, ONLY: initialize_molecule_type + USE input_types, ONLY: setup_parameters_type + USE integrator, ONLY: set_integrator + USE kinds, ONLY: int_size,& + wp => dp,& + wp_size => dp_size + USE md, ONLY: initialize_velocities,& + mdio_parameters_type,& + read_md_section,& + simulation_parameters_type,& + virial_type + USE md_environment_types, ONLY: init_md_env,& + md_environment_type,& + set_md_env,& + zero_virial + USE molecule_input, ONLY: read_molecule_section,& + read_setup_section,& + read_topology_section + USE molecule_types, ONLY: intra_parameters_type,& + molecule_structure_type,& + molecule_type,& + particle_node_type,& + topology_type + USE particle_types, ONLY: particle_prop_type,& + particle_type USE qs_build_ks_matrix, ONLY: qs_ks_env_type USE qs_energy, ONLY: qs_energies - USE qs_environment_types, ONLY: qs_environment_type + USE qs_environment_types, ONLY: get_qs_env,& + qs_environment_type USE qs_force, ONLY: qs_forces - USE termination, ONLY: stop_program + USE simulation_cell, ONLY: cell_type + USE structure_types, ONLY: structure_type,& + init_structure_type,& + set_structure_type + USE termination, ONLY: stop_memory,& + stop_program USE timings, ONLY: timeset,& timestop + USE unit, ONLY: unit_convert_type,& + set_units + USE velocity_verlet_control, ONLY: velocity_verlet IMPLICIT NONE @@ -68,7 +110,28 @@ CONTAINS ! *** Local variables *** - INTEGER :: handle + TYPE(extended_system_type) :: extended_type + TYPE(intra_parameters_type) :: intra_param + TYPE(md_environment_type) :: md_env + TYPE(mdio_parameters_type) :: mdio + TYPE(setup_parameters_type) :: setup + TYPE(topology_type) :: topo + TYPE(unit_convert_type) :: units + TYPE(virial_type) :: virial + + TYPE(cell_type), POINTER :: cell + TYPE(simulation_parameters_type), POINTER :: simpar + + TYPE(molecule_structure_type), DIMENSION(:), POINTER :: molecule + TYPE(molecule_type), DIMENSION(:), POINTER :: mol_setup + TYPE(particle_node_type), DIMENSION(:), POINTER :: pnode + TYPE(particle_prop_type), DIMENSION(:), POINTER :: pstat + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(structure_type), DIMENSION(:), POINTER :: struc + REAL(wp), POINTER :: constant + INTEGER, POINTER :: itimes + + INTEGER :: handle,istat,natom,nmolecule ! --------------------------------------------------------------------------- @@ -76,9 +139,97 @@ CONTAINS CALL timeset(routine_name,"I","",handle) -! CALL qs_forces(qs_env,globenv,ks_env) -! CALL qs_energies(qs_env,globenv,ks_env) - CALL stop_program(routine,"MD for QUICKSTEP not yet implemeted") + NULLIFY (molecule) + NULLIFY (mol_setup) + NULLIFY (pnode) + NULLIFY (pstat) + + ALLOCATE (simpar,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"simpar",0) + ALLOCATE (struc(1),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"struc",0) + + CALL read_cntl_section(setup,globenv) + + CALL read_setup_section(mol_setup,setup,globenv) + + CALL read_molecule_section(mol_setup,setup,globenv) + + CALL read_topology_section(topo,setup,globenv) + + CALL read_md_section(simpar,globenv,mdio) + simpar%program = globenv%program_name + + CALL set_particle_properties(setup,mol_setup,pstat) + + CALL set_units(setup%unit_type,units) + CALL convert(units=units,simpar=simpar) + + CALL get_qs_env(qs_env=qs_env,& + cell=cell,& + particle_set=particle_set) + + natom = SIZE(particle_set) + + ALLOCATE (pnode(natom),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"pnode",int_size*natom) + + nmolecule = SIZE(mol_setup) + + ALLOCATE (molecule(nmolecule),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"molecule",int_size*nmolecule) + + CALL initialize_molecule_type(mol_setup,intra_param,pnode,& + particle_set,molecule,globenv) + + ALLOCATE (extended_type%nhc_part(1),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"nhc_part",0) + + CALL initialize_nhc_part(cell,simpar,molecule,mol_setup,globenv,& + extended_type%nhc_part(1)) + + ALLOCATE (extended_type%nhc_baro(1),STAT=istat) + IF (istat /= 0 ) CALL stop_memory(routine,"nhc_baro",0) + + CALL initialize_nhc_baro(simpar,globenv,extended_type%nhc_baro(1)) + + CALL initialize_npt_type(simpar,globenv,extended_type%npt_info,cell) + + IF ((simpar%read_type == "POS").OR.& + (simpar%read_type == "INIT")) THEN + CALL initialize_velocities(simpar,particle_set,globenv) + END IF + + CALL set_force_env(qs_env_type=qs_env,& + ks_env_type=ks_env) + + CALL init_md_env(md_env) + + CALL set_md_env(md_env=md_env,& + simpar=simpar,& + struc=struc,& + nhc_part=extended_type%nhc_part,& + nhc_baro=extended_type%nhc_baro,& + npt=extended_type%npt_info,& + cell=cell) + + ALLOCATE (md_env%itimes,STAT=istat) + IF (istat /= 0) CALL stop_memory (routine,"itimes",int_size) + ALLOCATE (md_env%constant,STAT=istat) + IF (istat /= 0) CALL stop_memory (routine,"constant",wp_size) + + CALL set_integrator(globenv) + + itimes => md_env%itimes + itimes = 0 + + CALL qs_forces(qs_env,globenv,ks_env) + + DO itimes=1,md_env%simpar%nsteps +print*,"Step: :",itimes," of ",md_env%simpar%nsteps + CALL velocity_verlet(md_env) + + END DO CALL timestop(0.0_wp,handle) diff --git a/src/velocity_verlet_control.F b/src/velocity_verlet_control.F index 75b46a97c0..fc1effccbb 100644 --- a/src/velocity_verlet_control.F +++ b/src/velocity_verlet_control.F @@ -118,6 +118,15 @@ SUBROUTINE velocity_verlet ( md_env ) END SELECT END SELECT CASE ( QUICKSTEP ) + SELECT CASE (md_env%simpar%ensemble) + CASE ("NVE") + CALL nve(md_env) + CASE ("NVT") + CALL nvt(md_env) + CASE DEFAULT + CALL stop_program("velocity_verlet",& + "integrator not implemented for QUICKSTEP") + END SELECT CASE ( KIM_GORDON ) SELECT CASE ( md_env % simpar % ensemble ) CASE ( 'NVE')