comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 22:25:32 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

More MD stuff for QS

Browse source 
svn-origin-rev: 992
This commit is contained in:
Matthias Krack 2002-11-22 10:52:08 +00:00
parent dbf2f16008
commit 016ff5bc12
8 changed files with 380 additions and 39 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/fist_nonbond_force.F
View file

@ -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

3
src/force_control.F
View file

@ -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 )

15
src/particle_types.F
View file

@ -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

19
src/qs_energy.F
View file

@ -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)

202
src/qs_force.F
View file

@ -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

4
src/qs_main.F
View file

@ -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)

165
src/qs_md.F
View file

@ -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)

9
src/velocity_verlet_control.F
View file

@ -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')