Skip to content

Commit

Permalink
Add NHC-GLOBAL test + refactor
Browse files Browse the repository at this point in the history
  • Loading branch information
@danielhollas
danielhollas committed Mar 27, 2022
1 parent 4867e0c commit c925420
Show file tree
Hide file tree
Showing 16 changed files with 2,880 additions and 170 deletions.
3 changes: 1 addition & 2 deletions src/geom_analysis.F90
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,13 @@
module mod_analyze_geometry
use mod_const, only: DP, ANG, AUtoFS
use mod_array_size, only: NDISTMAX
use mod_general, only: sim_time, nwalk
implicit none
private
public :: ndist, nang, ndih, shiftdih
public :: dist1, dist2, ang1, ang2, ang3, dih1, dih2, dih3, dih4
public :: print_distances, print_angles, print_dihedrals

! TODO: Make all these allocatable, instead of using NDISTMAX
integer, parameter :: NDISTMAX = 50
integer :: ndist = 0, dist1(NDISTMAX), dist2(NDISTMAX)
integer :: nang = 0, ang1(NDISTMAX), ang2(NDISTMAX), ang3(NDISTMAX)
integer :: ndih = 0, dih1(NDISTMAX), dih2(NDISTMAX), dih3(NDISTMAX), dih4(NDISTMAX)
Expand Down
136 changes: 11 additions & 125 deletions src/init.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,9 +58,11 @@ subroutine init(dt)
real(DP) :: r0_morse = -1, d0_morse = -1, k_morse = -1
real(DP) :: k = -1, r0 = -1
real(DP) :: kx = -1, ky = -1, kz = -1
! Initial temperature (read from namelist nhcopt)
real(DP) :: temp0 = -1
real(DP) :: masses(MAXTYPES)
real(DP) :: rans(10)
integer :: ipom, iw, iat, natom_xyz, imol, iost
integer :: ipom, iw, iat, natom_xyz, iost
integer :: shiftdihed
! Number of OpenMP processes, read from ABIN input
! WARNING: We do NOT use OMP_NUM_THREADS environment variable!
Expand Down Expand Up @@ -407,10 +409,8 @@ subroutine init(dt)
! TODO: rename this function
call init_mass(amg, amt)

#if ( __GNUC__ == 4 && __GNUC_MINOR__ >= 6 ) || __GNUC__ > 4
allocate (natmolt(natom))
allocate (nshakemol(natom))
#endif
natmolt = 0
natmolt(1) = natom ! default for global NHC thermostat
nshakemol = 0
Expand Down Expand Up @@ -727,51 +727,17 @@ subroutine check_inputsanity()
end if

if (irest == 1 .and. chveloc /= '') then
! write(*,*)'ERROR: Input velocities are not compatible with irest=1.'
write (*, *) 'WARNING: Input velocities from file'//trim(chveloc)//' will be ignored!'
write (*, *) 'WARNING: Input velocities from file '//trim(chveloc)//' will be ignored!'
write (*, *) 'Velocities will be taken from restart file because irest=1.'
! write(*,*)chknow
! if(iknow.ne.1) error=1
end if

!-----Check, whether input variables don't exceeds array limits
if (ntraj > ntrajmax) then
write (*, *) 'Maximum number of trajectories is:'
write (*, *) ntrajmax
write (*, *) 'Adjust variable ntrajmax in modules.f90'
error = 1
end if
if (nstate > nstmax) then
write (*, *) 'Maximum number of states is:'
write (*, *) nstmax
write (*, *) 'Adjust variable nstmax in modules.f90'
error = 1
end if
if (nchain > maxchain) then
write (*, *) 'Maximum number of Nose-Hoover chains is:'
write (*, *) maxchain
write (*, *) 'Adjust variable maxchain in modules.f90'
error = 1
end if
if (ndist >= ndistmax) then
write (*, *) 'Maximum number of bonds for printing is:'
write (*, *) ndistmax
write (*, *) 'Adjust variable ndistmax in modules.f90'
error = 1
end if
if (nang >= ndistmax) then
write (*, *) 'Maximum number of angles (nang) for printing is:'
write (*, *) ndistmax
write (*, *) 'Adjust variable ndistmax in modules.f90'
error = 1
end if
if (ndih >= ndistmax) then
write (*, *) 'Maximum number of dihedral angles (ndih) for printing is:'
write (*, *) ndistmax
write (*, *) 'Adjust variable ndistmax in modules.f90'
error = 1
end if
!----------HERE we check for errors in input.

if (pot /= '_cp2k_') then
if (nproc > nwalk) then
write (*, *) 'ERROR: Nproc greater than nwalk. That does not make sense.'
Expand Down Expand Up @@ -803,7 +769,7 @@ subroutine check_inputsanity()
error = 1
end if
if (integ /= 'euler' .and. integ /= 'rk4' .and. integ /= 'butcher') then
write (*, *) 'integ must be "euler","rk4" or "butcher".'
write (*, *) 'integ must be "euler", "rk4" or "butcher".'
error = 1
end if
if (integ /= 'butcher') then
Expand Down Expand Up @@ -890,15 +856,13 @@ subroutine check_inputsanity()
end if

if (ipimd == 1 .and. inose <= 0) then
write (*, *) 'You have to use thermostat with PIMD!(inose>=0)'
write (*, *) 'You have to use thermostat with PIMD! (inose>=0)'
write (*, *) chknow
if (iknow /= 1) error = 1
end if
if (ipimd < 0 .or. ipimd > 3) then
if (ipimd /= 5) then
write (*, *) 'ipimd has to be 0,1,2 or 3.'
error = 1
end if
if (ipimd < 0 .or. ipimd > 5) then
write (*, *) 'Invalid ipimd value'
error = 1
end if
if (ipimd == 5 .and. pot == '_tera_' .and. ntriplet_lz > 0) then
write (*, *) 'ERROR: Landau-Zener with Singlet-Triplet transitions not implemented with TeraChem over MPI.'
Expand All @@ -916,13 +880,6 @@ subroutine check_inputsanity()
write (*, *) 'ERROR: inormalmodes has to be 0, 1 or 2!'
error = 1
end if
if (readnhc == 1 .and. initNHC == 1 .and. irest == 1) then
write (*, *) 'Warning: Conflicting keywords readnhc and initNHC set to 1.'
write (*, *) 'Momenta from restart.xyz will be used.'
end if
if (readnhc == 1 .and. irest == 0) then
write (*, *) 'Ignoring readnhc=1 since irest=0.'
end if
if (inac > 2 .or. inac < 0) then
write (*, *) 'Parameter "inac" must be 0,1 or 2.' !be very carefull if you change this!
error = 1
Expand Down Expand Up @@ -988,34 +945,6 @@ subroutine check_inputsanity()
write (*, *) 'nac_accu1 < nac_accu2'
write (*, *) 'I will compute NACME only with default accuracy:', nac_accu1
end if
if (imasst /= 0 .and. imasst /= 1) then
write (*, *) 'Input error: imasst must be 1 or zero.'
error = 1
end if
if (imasst == 0 .and. ipimd == 1) then
write (*, *) 'PIMD simulations must use massive thermostat ( imasst=1)! '
error = 1
end if
if (imasst == 0 .and. nmolt <= 0) then
write (*, *) 'Number of molecules coupled to separate NH chains not specified!Set nmolt > 0.'
error = 1
end if
if (nmolt > natom) then
write (*, *) 'Input error: nmolt > natom, which is not possible. Consult the manual.'
error = 1
end if
if (imasst == 0) then
do imol = 1, nmolt
if (natmolt(imol) <= 0) then
write (*, *) 'Number of atoms in molecules not specified! Set array natmolt properly.'
error = 1
end if
end do
end if
if (inose < 0 .and. inose > 3) then
write (*, *) 'inose has to be 0,1,2 or 3.'
error = 1
end if
if (istage == 1 .and. ipimd /= 1) then
write (*, *) 'The staging transformation is only meaningful for PIMD'
error = 1
Expand All @@ -1025,33 +954,14 @@ subroutine check_inputsanity()
error = 1
end if
if (istage == 0 .and. ipimd == 1 .and. inose /= 2 .and. inormalmodes == 0) then
write (*, *) 'PIMD should be done with staging or normal mode transformation! Exiting...'
write (*, *) 'PIMD should be done with staging or normal mode transformation!'
write (*, *) chknow
if (iknow /= 1) error = 1
end if
if (istage == 1 .and. inose == 2) then
write (*, *) 'The staging transformation is not compatible with GLE thermostat.'
error = 1
end if
if (nyosh /= 1 .and. nyosh /= 3 .and. nyosh /= 7) then
write (*, *) 'Variable nyosh(order of Suzuki-Yoshiga scheme) must be 1,3 or 7'
error = 1
end if
if (nyosh <= 1 .and. inose == 1) then
write (*, *) 'It is strongly reccommended to use Suzuki-Yoshida scheme when using Nose-Hoover thermostat (nyosh=3 or 7).'
write (*, *) iknow, error, chknow
if (iknow /= 1) error = 1
end if
if (nrespnose < 3 .and. inose == 1) then
write (*, *) 'Variable nrespnose < 3! Assuming this is an error in input and exiting.'
write (*, *) 'Such low value would probably not produce stable results.'
write (*, *) chknow
if (iknow /= 1) error = 1
end if
if (nrespnose <= 0) then
write (*, *) 'Variable nrespnose must be positive integer'
error = 1
end if
if (irest /= 1 .and. irest /= 0) then
write (*, *) 'ERROR:irest has to be 1 or 0'
error = 1
Expand All @@ -1065,10 +975,6 @@ subroutine check_inputsanity()
write (*, *) 'Set imasst=1 and nmolt, natmolt and nshakemol accordingly.'
error = 1
end if
if (pot == 'mm' .and. iqmmm > 0) then
write (*, *) 'Pot="mm"is not compatible with iqmmm>0!'
error = 1
end if
if (iqmmm > 1) then
write (*, *) 'WARNING: QMMM is higly experimental at this point. Use with care!'
write (*, *) chknow
Expand All @@ -1079,26 +985,6 @@ subroutine check_inputsanity()
error = 1
end if

if (inose == 1 .and. imasst == 0) then
ipom = 0
do iat = 1, nmolt
ipom = ipom + natmolt(iat)
end do
if (ipom /= natom) then
write (*, *) 'Number of atoms in thermostated molecules(natmolt) doesnt match with natom.'
write (*, *) 'This is probably mistake in input.Exiting...'
write (*, *) chknow
if (iknow /= 1) error = 1
end if
end if

if (temp < 1 .and. inose >= 1) then
write (*, *) 'WARNING!:Temperature below 1K. Are you sure?'
write (*, *) 'This is probably mistake in input.Exiting...'
write (*, *) chknow
if (iknow /= 1) error = 1
end if

if (iremd == 1) then
write (chout, '(A,I2.2)') 'movie.xyz.', my_rank
else
Expand Down
3 changes: 1 addition & 2 deletions src/modules.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,7 @@ module mod_array_size
use mod_const, only: DP
implicit none
public
integer, parameter :: MAXCHAIN = 10, MAXTYPES = 10
integer, parameter :: NBINMAX = 2000, NDISTMAX = 30
integer, parameter :: MAXTYPES = 10
integer, parameter :: NSTMAX = 50, NTRAJMAX = 1
save
end module mod_array_size
Expand Down
Loading

0 comments on commit c925420

Please sign in to comment.