#!/usr/bin/env bash # These files will be patched and must exist files=("./doc/manual/manual.tex" "./scripts/configure.sh" "./src/anneal_utils.mod.F90" "./src/control_def_utils.mod.F90" "./src/control_pri_utils.mod.F90" "./src/control_utils.mod.F90" "./src/cotr.mod.F90" "./src/cpmd.F90" "./src/detdof_utils.mod.F90" "./src/eextern_utils.mod.F90" "./src/finalp_utils.mod.F90" "./src/forcedr_driver.mod.F90" "./src/forces_diag_utils.mod.F90" "./src/forces_driver.mod.F90" "./src/geofile_utils.mod.F90" "./src/initrun_driver.mod.F90" "./src/loadpa_utils.mod.F90" "./src/md_driver.mod.F90" "./src/mdmain_utils.mod.F90" "./src/mm_init_utils.mod.F90" "./src/mp_interface.mod.F90" "./src/mts_utils.mod.F90" "./src/nosalloc_utils.mod.F90" "./src/nosepa_utils.mod.F90" "./src/pbc_utils.mod.F90" "./src/pcgrad_driver.mod.F90" "./src/potfor_utils.mod.F90" "./src/printp_utils.mod.F90" "./src/proppt_utils.mod.F90" "./src/ratom_utils.mod.F90" "./src/rv30_utils.mod.F90" "./src/rwfopt_utils.mod.F90" "./src/setirec_utils.mod.F90" "./src/setsys_utils.mod.F90" "./src/shake_utils.mod.F90" "./src/SOURCES" "./src/sysin_utils.mod.F90" "./src/system.mod.F90" "./src/updrho_utils.mod.F90" "./src/wrener_utils.mod.F90" "./src/wrgeo_utils.mod.F90" "./src/wv30_utils.mod.F90") # Check if files exist for file in ${files[@]} do if [ ! -w "$file" ] then echo "ERROR: $file does not exist or is not writable. Are we in the CPMD root directory?" exit fi done # These files will be created and must not exist new_files=("./configure/MIMIC-LINUX-GCC" "./configure/MIMIC-LINUX-INTEL" "./src/constraint.mod.F90" "./src/constraint_utils.mod.F90" "./src/mimic_wrapper.mod.F90") # Check that files do not exist for file in ${new_files[@]} do if [ -e "$file" ] then echo "ERROR: $file already exists. Has CPMD already been patched?" exit fi done # Create patches cat << EOF > patch.00 diff ./configure/MIMIC-LINUX-GCC ./configure/MIMIC-LINUX-GCC 0a IRAT=2 FC='mpif90' CC='mpicc' LD='mpif90' CPP='cpp -P -traditional' CPPFLAGS=' -D__Linux -D__HAS_FFT_FFTW3 -D__PARALLEL -D__HAS_SIZEOF ' if [ \$debug ]; then FFLAGS=' -g -pg -Og -ffree-line-length-none ' CFLAGS=' -g -pg -Og ' else FFLAGS=' -O3 -march=native -mtune=native ' CFLAGS=' -O3 -march=native -mtune=native ' fi if [ \$omp ]; then FFLAGS=\${FFLAGS}' -fopenmp ' OMP3_DISABLED='false' LIBS=' -L/usr/lib64 -llapack -lopenblaso -lfftw3_omp -lfftw3 -lm ' else LIBS=' -L/usr/lib64 -llapack -lopenblas -lfftw3 -lm ' fi LFLAGS=\${LIBS} FFLAGS=\${FFLAGS}" -I/usr/include -fallow-argument-mismatch -ffree-line-length-none" LFLAGS=\${LFLAGS}" -L/usr/lib64 " . EOF cat << EOF > patch.01 diff ./configure/MIMIC-LINUX-INTEL ./configure/MIMIC-LINUX-INTEL 0a IRAT=2 FC='mpiifort' CC='mpiicc' LD='mpiifort' CPP='fpp' CPPFLAGS=' -D__Linux -D__HAS_FFT_FFTW3 -D__PARALLEL -D__HAS_SIZEOF -D__HAS_BF_STREAM_IO ' AR='xiar -r' RANLIB='ranlib' if [ \$debug ]; then FFLAGS=' -g -O0 -check all -warn all -traceback ' CFLAGS=' -g -O0 -check all -warn all -traceback ' else FFLAGS=' -xHost -O3 ' CFLAGS=' -xHost -O3 ' fi if [ \$omp ]; then FFLAGS=\${FFLAGS}' -qopenmp ' OMP3_DISABLED='false' LIBS=' -qmkl=parallel -lm ' else LIBS=' -qmkl=sequential -lm ' fi LFLAGS=\${LIBS} FFLAGS=\${FFLAGS}" -I/usr/include " LFLAGS=\${LFLAGS}" -L/usr/lib64 " . EOF cat << EOF > patch.02 diff ./doc/manual/manual.tex ./doc/manual/manual.tex 11986a \bibitem{Olsen19} J.~M.~H.~Olsen, V.~Bolnykh, S.~Meloni, E.~Ippoliti, M.~P.~Bircher, P.~Carloni, U.~Rothlisberger, J. Chem. Theory Comput. {\bf 15}, 3810--3823 (2019). DOI: \href{https://doi.org/10.1021/acs.jctc.9b00093}{10.1021/acs.jctc.9b00093} \bibitem{Bolnykh19} V.~Bolnykh, J.~M.~H.~Olsen, S.~Meloni, E.~Ippoliti, M.~P.~Bircher, P.~Carloni, U.~Rothlisberger, J. Chem. Theory Comput. {\bf 15}, 5601--5613 (2019). DOI: \href{https://doi.org/10.1021/acs.jctc.9b00424}{10.1021/acs.jctc.9b00424} \bibitem{Weinback05} Y.~Weinback, R.~Elber, J. Comput. Phys. {\bf 209}, 193--206 (2005). DOI: \href{https://www.sciencedirect.com/science/article/pii/S0021999105001828}{10.1016/j.jcp.2005.03.015} . 10391c PROGRAMER} (during startup of WAVEFUNCTION or GEOMETRY optimization)} . 5455a \keyword{TOTAL SCF}{}{}{}{\&MIMIC} \desc{Force MiMiC to request MM energies from client programs in the beginning of an SCF loop (it is always retrieved after each completed SCF). Without this keyword, the energy during an SCF loop will be missing the MM energy.\\\\ {\bf Note} that using this keyword will trigger additional communication with the clients and will also prevent simultaneous calculation of QM and MM contributions.} . 5397c value and two more numbers are read, the ramping target temperature in Kelvin and the . 5222a \keyword{SUBSYSTEM FACTORS}{}{}{}{\&MIMIC} \desc{Define the factors multiplying the energy, forces, etc. of each subsystem. Useful only for certain types of simulations, e.g., subtractive QM/MM.\\\\ The number of factors, which must be equal to the number of clients, is read from the next line. Then the factors are read from the subsequent lines.\\\\ {\bf Default} is to use 1.0 for all subsystems.} . 5054a \keyword{SHAKE_\CG\_ITER}{}{\&CPMD} \desc{The maximum number of iterations of the minimizer is read from the next line (see \refkeyword{NEW CONSTRAINTS}).} \keyword{SHAKE\_MAXSTEP}{}{\&CPMD} \desc{The maximum number of shake iterations is read from the next line (see \refkeyword{NEW CONSTRAINTS}).} . 4276a \keyword{PATHS}{}{}{}{\&MIMIC} \desc{Define the paths to the working directories from which the client programs are launched. The number of clients (paths) is read from the next line followed by the paths.} . 4245a \keyword{OVERLAPS}{}{}{}{\&MIMIC} \desc{Define the mapping of the atoms that are present in more than one subsystem. For QM/MM this will typically be the atoms in the QM subsystem which are used by both the QM and MM client codes.\\\\ The format is the following: \begin{center} \begin{tabular}{llll} \multicolumn{4}{l}{num\_overlap\_pairs} \\\\ client\_id & client\_atom\_id & host\_id & host\_atom\_id \\\\ \end{tabular} \end{center} where num\_overlap\_pairs is the number of pairs of atom overlaps (and the number of succeeding lines). Each of the following lines define the mapping of atom pairs. Here client\_id is the id of the client code, which corresponds to the number of the line in which this client appears in the \refkeyword{PATHS} section plus one (e.g., the first line in \refkeyword{PATHS} has client\_id = 2, the second line has client\_id = 3, etc.), client\_atom\_id is the id of the atom in the given client, host\_id is the id of the MD driver (currently CPMD and its value is always 1), and host\_atom\_id is the id of the atom in the MD driver (which follows the ordering of the atoms as provided in the \&ATOMS \ldots\ \&END section).\\\\ An example input could be:\\\\ OVERLAPS\\\\ 3\\\\ 2 1 1 2\\\\ 2 2 1 3\\\\ 2 3 1 1\\\\ which means that the atom with id 1 in the client code with id 2 overlaps with atom 2 in the host code (CPMD), the atom with id 2 in the client code with id 2 overlaps with atom 3 in the host code, and the atom with id 3 in the client code with id 2 overlaps with atom 1 in the host code.} . 4117c the GROMOS ordering in case of classical or QM/MM runs. In the case of MiMiC-based QM/MM, only CPMD ordering can be used. Multiple . 3966a \keyword{NEW CONSTRAINTS}{}{[PBICGSTAB]}{}{\&CPMD} \desc{This keyword activates the alternative implementation of the shake algorithm described in~\cite{Weinback05}. This implementation does not involve matrix inversion operations and is therefore (much) faster than the default one. It can be used in pure QM simulations and in \refkeyword{MiMiC}-based QM/MM simulations. At each iteration of the shake algorithm, a minimization is performed using a preconditioned conjugate gradient (PCG) method, unless the PBICGSTAB optional keyword is specified, in which case, a stabilized precondictioned bi-conjugate gradient (PBiCGStab) method is adopted. In general, we recommend using the default PCG minimizer, which is slightly faster. The \refkeyword{SHAKE\_MAXSTEP} and \refkeyword{SHAKE\_CG\_ITER} keywords can be used to set the maximum number of shake iterations and the maximum number of steps of the PCG (or PBiCGStab) minimizer, respectively.\\\\ {\bf Note:} Currently only a serial version of the algorithm is implemented.} . 3942a \keyword{MULTIPOLE ORDER}{}{}{}{\&MIMIC} \desc{Define the order of the multipole expansion used in the \refkeyword{LONG-RANGE COUPLING} where electrostatic QM/MM interactions are split into short- and long-range parts.\\\\ {\bf Default} is to use a second-order multipole expansion.} . 3841a \keyword{MIMIC}{}{}{}{\&CPMD} \desc{Activate multiscale simulations using MiMiC. See \url{https://mimic-project.org/} for more information on how to setup and run MiMiC-based simulations and Refs.~\cite{Olsen19} and \cite{Bolnykh19} for more information about MiMiC.} . 3742a \keyword{LONG-RANGE COUPLING}{}{}{}{\&MIMIC} \desc{Activate the long-range coupling that partitions the electrostatic QM/MM interactions into short- and long-range parts. The short-range interactions are computed using a damped Coulomb potential whereas a multipole expansion of the QM electrostatic potential is used to treat the long-range interactions.\\\\ The cutoff distance for the partitioning is controlled by the \refkeyword{CUTOFF DISTANCE} keyword, the order of the multipole expansion is controlled by the \refkeyword{MULTIPOLE ORDER} keyword, and the method and frequency of sorting atoms into short- and long-range parts is controlled by the \refkeyword{FRAGMENT SORTING} keyword.} . 2980a \keyword{FRAGMENT SORTING}{\{CENTROID, CENTER-OF-MASS, CENTER-OF-CHARGE, MINIMUM DISTANCE, ATOM-WISE\}}{}{[UPDATE]}{\&MIMIC} \desc{Define the method used to sort atoms into short- and long-range domains. It is used only when \refkeyword{LONG-RANGE COUPLING} is active.\\\\ The CENTROID, CENTER-OF-MASS, CENTER-OF-CHARGE, and MINIMUM DISTANCE methods use the distance between centroids, centers of mass, centers of charge, and minimum interatomic distances, respectively, between the QM subsystem and fragments in the MM subsystem. The ATOM-WISE method uses the distance between the centroid of the QM subsystem and the individual atoms in the environment. The frequency of the sorting updates can be set using the {\bf UPDATE} keyword. When present the frequency (in steps) is read from the next line.\\\\ {\bf Default} is to use CENTROID and to update the sorting in each step.} . 2608a \keyword{DISABLE CONSTRAINTS}{}{}{}{\&MIMIC} \desc{This keyword instructs CPMD to ignore constraints defined by client programs in MIMIC-based multiscale simulations. Only the constraints defined in the CPMD input file will be enforced.} . 2425a \keyword{CUTOFF DISTANCE}{}{}{}{\&MIMIC} \desc{Set the cutoff distance for the partitioning of the MM subsystem into short- and long-range domains. The distance (in bohr) is read from the next line. It is used for the \refkeyword{LONG-RANGE COUPLING} in MiMiC-based QM/MM simulations.\\\\ {\bf Default} is to use a large cutoff distance placing all atoms in the short-range domain.} . 1930a \keyword{BOX}{}{}{}{\&MIMIC} \desc{Define the dimensions of the outer simulation box in MiMiC-based multiscale simulations. Values are read from the next line.\\\\ {\bf Note:} currently only orthorombic simulation boxes are supported.} . 1762c time step. The scaling factor is read from the next line.\\\\ The additional DUAL keyword can be used in combination with the IONS keyword in \refkeyword{MIMIC}-based QM/MM simulations. It enables defining two scaling factors, the first that will be applied to the QM atoms only and the second that will be applied to the MM atoms only. The two scaling factors are read from the next line.} . 1760c \keyword{ANNEALING}{\{IONS,ELECTRONS,CELL\}}{DUAL}{}{\&CPMD} . 1662a \subsubsection[\&MIMIC ... \&END]{\&MIMIC \$\ldots\$ \&END} % \options{}{}{} \refkeyword{BOX}\options{}{}{} \refkeyword{CUTOFF DISTANCE}\options{}{}{} \refkeyword{DISABLE CONSTRAINTS}\options{}{}{} \refkeyword{FRAGMENT SORTING}\options{\{CENTROID, CENTER-OF-MASS, \\\\ \phantom{FRAGMENT SORTING \{} CENTER-OF-CHARGE, MINIMUM DISTANCE, \\\\ \phantom{FRAGMENT SORTING \{} ATOM-WISE\}}{[UPDATE]}{} \refkeyword{LONG-RANGE COUPLING}\options{}{}{} \refkeyword{MULTIPOLE ORDER}\options{}{}{} \refkeyword{OVERLAPS}\options{}{}{} \refkeyword{PATHS}\options{}{}{} \refkeyword{SUBSYSTEM FACTORS}\options{}{}{} \refkeyword{TOTAL SCF}\options{}{}{} % . 1255a \refkeyword{SHAKE\_CG\_ITER}\options{}{}{} \refkeyword{SHAKE\_MAXSTEP}\options{}{}{} . 1216a \refkeyword{NEW CONSTRAINTS}\options{}{[PBICGSTAB]}{} . 1209a \refkeyword{MIMIC}\options{}{}{} . 1167a \refkeyword{DISABLE CONSTRAINTS}\options{}{}{} . 1136c \refkeyword{ANNEALING}\options{\{IONS,ELECTRONS,CELL\}}{[DUAL]}{} . 1119a \\\\ \> \&MIMIC ... \> \&END \> \$\leftrightarrow\$ \> Parameters for multiscale simulations using MiMiC. \\\\ \> This section is only evaluated if the \refkeyword{MIMIC} keyword is given in the \&CPMD section.\\\\ . 659a Coupling to the MiMiC multiscale modeling framework (see \url{https://mimic-project.org/}) is added which currently enables QM/MM simulations using GROMACS as the MM engine. . EOF cat << EOF > patch.03 diff ./scripts/configure.sh ./scripts/configure.sh 854a SkipInclude["mimic_precision"] = 0; SkipInclude["mimic_constants"] = 0; SkipInclude["mimic_types"] = 0; SkipInclude["mimic_tensors"] = 0; SkipInclude["mimic_main"] = 0; . 243a if [ \$mimic ]; then mkdir mimic_test cd mimic_test echo "program test" > test.f90 echo "use mimic_main" >> test.f90 echo "end program test" >> test.f90 \${FC} test.f90 \${FFLAGS} \${LFLAGS} -lmimic -lmimiccommf -lmimiccomm if [ \$? != 0 ] then cd .. rm -r mimic_test echo "MiMiC was not found on your system!" >&2 exit 1 fi cd .. rm -r mimic_test CPPFLAGS=\${CPPFLAGS}' -D__MIMIC' LFLAGS=\${LFLAGS}' -lmimic -lmimiccommf -lmimiccomm -lstdc++ ' fi . 132a -mimic|-m) mimic=1 echo "** Enabling MiMiC interface (if MiMiC is available)" >&2 ;; . 60c -mimic Enable interface to MiMiC (requires MPI and that MiMiC is available) -qmmm Generates a makefile for QMMM . EOF cat << EOF > patch.04 diff ./src/anneal_utils.mod.F90 ./src/anneal_utils.mod.F90 60,66c DO i = 1, mimic_control%num_quantum_atoms ia = iatpt(1,i) is = iatpt(2,i) velp(1,ia,is) = alfap * velp(1,ia,is) velp(2,ia,is) = alfap * velp(2,ia,is) velp(3,ia,is) = alfap * velp(3,ia,is) ENDDO alfap = cntr%anneal_factors(2)**(0.25_real_8) #if defined(__VECTOR) !\$omp parallel do private(I,IS,IA) #else !\$omp parallel do private(I,IS,IA) schedule(static) #endif DO i = mimic_control%num_quantum_atoms + 1, mimic_control%num_atoms ia = iatpt(1,i) is = iatpt(2,i) velp(1,ia,is) = alfap * velp(1,ia,is) velp(2,ia,is) = alfap * velp(2,ia,is) velp(3,ia,is) = alfap * velp(3,ia,is) ENDDO ELSE alfap=cntr%anneri**(0.25_real_8) #if defined(__VECTOR) !\$omp parallel do private(I,IS,IA) #else !\$omp parallel do private(I,IS,IA) schedule(static) #endif DO i=1,ions1%nat ia=iatpt(1,i) is=iatpt(2,i) velp(1,ia,is)=alfap*velp(1,ia,is) velp(2,ia,is)=alfap*velp(2,ia,is) velp(3,ia,is)=alfap*velp(3,ia,is) ENDDO END IF . 58c !\$omp parallel do private(I,IS,IA) schedule(static) . 56c !\$omp parallel do private(I,IS,IA) . 54c IF (cntl%anneal_dual) THEN alfap = cntr%anneal_factors(1)**(0.25_real_8) . 7a USE mimic_wrapper, ONLY: mimic_control . EOF cat << EOF > patch.05 diff ./src/constraint.mod.F90 ./src/constraint.mod.F90 0a module constraint USE kinds, ONLY: real_8 implicit none !> Type used to store sparse constraint matrix !> @author V. Bolnykh - RWTH Aachen/Cyprus Institute type constraint_matrix ! Values stored in a sparse matrix real(real_8), dimension(:), allocatable :: vals ! Number of non-zero elements of the current row integer, dimension(:), allocatable :: nz_row ! Index of a non-zero element in the row integer, dimension(:), allocatable :: r_index ! IDs of atoms in a constraint integer, dimension(:, :), allocatable :: atom_id ! Species and atom IDs in the constraint (needed for TAU ! mapping) integer, dimension(:, :, :), allocatable :: tau_map ! Reference value for a constraint real(real_8), dimension(:), allocatable :: ref contains procedure :: mat_vec generic :: operator(*) => mat_vec end type constraint_matrix !> Interface of the constraint_matrix constructor interface constraint_matrix module procedure new_matrix end interface contains !> Matrix-vector multiplication function mat_vec(m, v) result(r) !> input matrix class(constraint_matrix), intent(in) :: m !> input vector real(real_8), dimension(:), intent(in) :: v !> resulting vector real(real_8), dimension(:), allocatable :: r integer :: i, j allocate(r(size(v))) r = 0.0_real_8 !\$OMP PARALLEL DO PRIVATE(i, j) do i = 1, size(m%nz_row) - 1 do j = m%nz_row(i - 1) + 1, m%nz_row(i) r(i) = r(i) + m%vals(j) * v(m%r_index(j)) end do end do !\$OMP END PARALLEL DO end function !> Constructor of the sparse matrix function new_matrix(vals, nz_row, r_index, atom_id, tau_map, ref) real(real_8), dimension(:), intent(in) :: vals integer, dimension(0:), intent(in) :: nz_row integer, dimension(:), intent(in) :: r_index integer, dimension(:,:), intent(in) :: atom_id integer, dimension(:,:,:), intent(in) :: tau_map real(real_8), dimension(:), intent(in) :: ref type(constraint_matrix) :: new_matrix integer :: i, j allocate(new_matrix%vals(size(vals))) allocate(new_matrix%nz_row(0:size(nz_row) - 1)) allocate(new_matrix%r_index(size(r_index))) allocate(new_matrix%atom_id(size(atom_id, 1), size(atom_id, 2))) allocate(new_matrix%tau_map(size(tau_map, 1), size(tau_map, 2), & size(tau_map, 3))) allocate(new_matrix%ref(size(ref))) !\$OMP PARALLEL PRIVATE(i, j) !\$OMP DO do i = 1, size(vals) new_matrix%vals(i) = vals(i) end do !\$OMP END DO !\$OMP DO do i = 0, size(nz_row) - 1 new_matrix%nz_row(i) = nz_row(i) end do !\$OMP END DO !\$OMP DO do i = 1, size(r_index) new_matrix%r_index(i) = r_index(i) end do !\$OMP END DO !\$OMP DO do i = 1, size(atom_id, 2) new_matrix%atom_id(:, i) = atom_id(:, i) end do !\$OMP END DO !\$OMP DO do i = 1, size(tau_map, 3) do j = 1, size(tau_map, 2) new_matrix%tau_map(:, j, i) = tau_map(:, j, i) end do end do !\$OMP END DO !\$OMP DO do i = 1, size(ref) new_matrix%ref(i) = ref(i) end do !\$OMP END DO !\$OMP END PARALLEL end function end module constraint . EOF cat << EOF > patch.06 diff ./src/constraint_utils.mod.F90 ./src/constraint_utils.mod.F90 0a module constraint_utils USE system, ONLY: cntl, cnti USE parac, ONLY: paral USE machine, ONLY: m_walltime USE constraint USE constr_utils, ONLY: funcr, diffr USE dum2_utils, ONLY: dum2 USE error_handling, ONLY: stopgm USE fillc_utils, ONLY: fillc USE kinds, ONLY: real_8 USE timer, ONLY: tihalt,& tiset implicit none ! stretch constraint type integer, parameter :: TYPE_STRETCH = 1 contains !> Initialize constraints matrix function build_constraints(ntcnst, cval, na, nsp) result(matrix) ! CPMD constraint map integer, dimension(:,:), intent(in) :: ntcnst ! Reference values for constraints real(real_8), dimension(:), intent(in) :: cval ! Number of atoms per species integer, dimension(:), intent(in) :: na ! Number of atomic species integer, intent(in) :: nsp ! Resulting sparse matrix type(constraint_matrix) :: matrix integer :: i, j, iat, isp integer :: nconstr, n_coupl integer :: ctype, ctype2, ia, ia2, ib, ib2, id real(real_8), dimension(:), allocatable :: vals integer, dimension(:), allocatable :: nz_row integer, dimension(:), allocatable :: nz_column, temp_nz integer, dimension(:,:), allocatable :: ids integer, dimension(:,:,:), allocatable :: tau_map integer :: accumulator ! Count constraints nconstr = 0 do i = 1, size(ntcnst, 2) ctype = ntcnst(1, i) if (ctype == TYPE_STRETCH) then nconstr = nconstr + 1 else CALL stopgm('BUILD_CONSTRAINTS','UNSUPPORTED CONSTRAINT TYPE',& __LINE__,__FILE__) end if end do ! i allocate(nz_row(0:nconstr)) !FIXME update for other types of constraints? allocate(ids(2, nconstr)) allocate(tau_map(2, 2, nconstr)) nz_row = 0 ids = 0 tau_map = 0 ! Check coupling between constraints n_coupl = 0 accumulator = 0 do i = 1, size(ntcnst, 2) ctype = ntcnst(1, i) ia = ntcnst(2, i) ib = ntcnst(3, i) ids(1, i) = ia ids(2, i) = ib do j = 1, size(ntcnst, 2) ctype2 = ntcnst(1, j) if (ctype2 == TYPE_STRETCH) then ia2 = ntcnst(2, j) ib2 = ntcnst(3, j) if (ia == ia2 .or. ib == ib2 .or. & ia == ib2 .or. ib == ia2) then accumulator = accumulator + 1 n_coupl = n_coupl + 1 allocate(temp_nz(n_coupl)) temp_nz = 0 if (allocated(nz_column)) then temp_nz(1:n_coupl - 1) = nz_column deallocate(nz_column) end if call move_alloc(temp_nz, nz_column) nz_column(n_coupl) = j end if ! overlap end if ! ctype end do ! j id = 0 do isp = 1, nsp do iat = 1, na(isp) id = id + 1 if (id == ia) then tau_map(1, 1, i) = isp tau_map(2, 1, i) = iat end if if (id == ib) then tau_map(1, 2, i) = isp tau_map(2, 2, i) = iat end if end do end do nz_row(i) = accumulator end do ! i allocate(vals(size(nz_column))) vals = 0.0_real_8 matrix = new_matrix(vals, nz_row, nz_column, ids, tau_map, cval) end function build_constraints !> Update constraint values and their gradients subroutine update_constraints(matrix, grad, diff, tau) !> Constraint matrix - needed to determine connectivity class(constraint_matrix), intent(inout) :: matrix !> Constraint gradients, output real(real_8), dimension(:, :), intent(out) :: grad !> Constraint violation, output real(real_8), dimension(:), intent(out) :: diff !> Coordinate matrix real(real_8), dimension(:,:,:), intent(in) :: tau character(*), PARAMETER :: procedureN = 'update_constraints' integer :: isub integer :: i, j integer :: ia, ib integer :: cid real(real_8) :: dist, dist_sq real(real_8), dimension(3) :: r_a, r_b real(real_8), dimension(6) :: d_sigma call tiset(procedureN,isub) grad = 0.0_real_8 !\$OMP PARALLEL DO PRIVATE(i, ia, ib, r_a, r_b, dist_sq, d_sigma) do i = 1, size(matrix%nz_row) - 1 ia = matrix%atom_id(1, i) ib = matrix%atom_id(2, i) call fillc(ia, tau, r_a) call fillc(ib, tau, r_b) call funcr(diff(i), dist_sq, matrix%ref(i), r_a, r_b) call diffr(d_sigma, r_a, r_b) grad(:, i) = grad(:, i) + d_sigma end do ! i !\$OMP END PARALLEL DO call tihalt(procedureN,isub) end subroutine update_constraints !> Fill symmetric constraint matrix subroutine fill_matrix(matrix, grad, dt, dtm) !> Matrix to fill (SHOULD BE INITIALIZED) class(constraint_matrix), intent(inout) :: matrix !> Constraint gradients real(real_8), dimension(:,:), intent(in) :: grad !> Ionic timestep real(real_8), intent(in) :: dt !> dt / (2 * mi) real(real_8), dimension(:), intent(in) :: dtm character(*), PARAMETER :: procedureN = 'fill_matrix' integer :: isub integer :: i, j, k, is, ia, pos1, pos2, at1, at2, cid real(real_8), dimension(:), allocatable :: vals call tiset(procedureN,isub) allocate(vals(size(matrix%vals))) vals = 0.0_real_8 !\$OMP PARALLEL DO PRIVATE(i, j, cid, at1, pos1, at2, pos2) REDUCTION(+:vals) do i = 1, size(matrix%nz_row) - 1 do j = matrix%nz_row(i - 1) + 1, matrix%nz_row(i) cid = matrix%r_index(j) do at1 = 1, size(matrix%atom_id, 1) pos1 = (at1 - 1) * 3 + 1 do at2 = 1, size(matrix%atom_id, 1) pos2 = (at2 - 1) * 3 + 1 if (matrix%atom_id(at1, i) == matrix%atom_id(at2, cid)) then vals(j) = vals(j) + & dt * dtm(matrix%tau_map(1, at1, i)) * & dot_product(grad(pos1:pos1+2, i), grad(pos2:pos2+2, cid)) end if end do end do end do end do !\$OMP END PARALLEL DO deallocate(matrix%vals) call move_alloc(vals, matrix%vals) call tihalt(procedureN,isub) end subroutine fill_matrix !> Fill asymmetric constraint matrix subroutine fill_matrix_new(matrix, gradl, gradr, dt, dtm) !> Matrix to fill (SHOULD BE INITIALIZED) class(constraint_matrix), intent(inout) :: matrix !> Constraint gradients real(real_8), dimension(:,:), intent(in) :: gradl, gradr !> Ionic timestep real(real_8), intent(in) :: dt !> dt / (2 * mi) real(real_8), dimension(:), intent(in) :: dtm character(*), PARAMETER :: procedureN = 'fill_matrix_new' integer :: isub integer :: i, j, k, is, ia, pos1, pos2, at1, at2, cid real(real_8), dimension(:), allocatable :: vals call tiset(procedureN,isub) allocate(vals(size(matrix%vals))) vals = 0.0_real_8 !\$OMP PARALLEL DO PRIVATE(i, j, cid, at1, pos1, at2, pos2) REDUCTION(+:vals) do i = 1, size(matrix%nz_row) - 1 do j = matrix%nz_row(i - 1) + 1, matrix%nz_row(i) cid = matrix%r_index(j) do at1 = 1, size(matrix%atom_id, 1) pos1 = (at1 - 1) * 3 + 1 do at2 = 1, size(matrix%atom_id, 1) pos2 = (at2 - 1) * 3 + 1 if (matrix%atom_id(at1, i) == matrix%atom_id(at2, cid)) then vals(j) = vals(j) + & dt * dtm(matrix%tau_map(1, at1, i)) * & dot_product(gradl(pos1:pos1+2, i), gradr(pos2:pos2+2, cid)) end if end do end do end do end do !\$OMP END PARALLEL DO deallocate(matrix%vals) call move_alloc(vals, matrix%vals) call tihalt(procedureN,isub) end subroutine fill_matrix_new !> Perform SHAKE constraint solver subroutine new_shake(matrix, xlagr, tau0, taup, velp, dt, dtm) !> Constraint matrix class(constraint_matrix), intent(inout) :: matrix !> Lagrange multiplier array real(real_8), dimension(:) :: xlagr !> CPMD coordinate matrix real(real_8), dimension(:,:,:) :: tau0 real(real_8), dimension(:,:,:) :: taup !> CPMD velocity matrix real(real_8), dimension(:,:,:) :: velp !> Ionic timestep real(real_8), intent(in) :: dt !> dt / (2 * mi) real(real_8), dimension(:), intent(in) :: dtm integer :: i, at, atid, spid, cat, pos real(real_8), dimension(:), allocatable :: sigma real(real_8), dimension(:, :), allocatable :: grad_sigma real(real_8), dimension(:, :), allocatable :: grad_current real(real_8), dimension(:,:,:), allocatable :: tscr real(real_8), parameter :: threshold = 1.0e-7 real(real_8) :: time1, time2 real(real_8), external :: dnrm2 character(*), PARAMETER :: procedureN = 'new_shake' integer :: isub call tiset(procedureN,isub) time1 = m_walltime() ! Allocate temporary storage allocate(sigma(size(matrix%nz_row)-1)) allocate(grad_sigma(6, size(matrix%nz_row)-1)) allocate(grad_current(6, size(matrix%nz_row)-1)) allocate(tscr(size(tau0, 1), size(tau0, 2), size(tau0, 3))) sigma = 0.0_real_8 grad_sigma = 0.0_real_8 grad_current = 0.0_real_8 tscr = 0.0_real_8 call dum2(tau0, tscr) call update_constraints(matrix, grad_sigma, sigma, tscr) do i = 1, cnti%shake_maxstep call dum2(taup,tscr) call update_constraints(matrix, grad_current, sigma, tscr) call fill_matrix(matrix, grad_current, dt, dtm) if (dnrm2(size(sigma), sigma, 1) < threshold) then time2 = m_walltime() if (paral%io_parent) then write(6,'(8x,a,i3,a,t58,f8.2)') "SHAKE CONVERGED AFTER ", i, & " ITERATIONS. TCPU = ", & (time2 - time1) * 0.001_real_8 end if exit end if if (cntl%pbicgstab) then call pbicgstab_solve(matrix, xlagr, sigma, threshold) else call pcg_solve(matrix, xlagr, sigma, threshold) endif do at = 1, size(matrix%nz_row) - 1 do cat = 1, size(matrix%atom_id, 1) spid = matrix%tau_map(1, cat, at) atid = matrix%tau_map(2, cat, at) pos = (cat - 1) * 3 + 1 taup(:, atid, spid) = taup(:, atid, spid) - dtm(spid) * dt * xlagr(at) * grad_sigma(pos : pos + 2, at) velp(:, atid, spid) = velp(:, atid, spid) - dtm(spid) * xlagr(at) * grad_sigma(pos : pos + 2, at) end do end do if (i == cnti%shake_maxstep) then if (paral%io_parent) then write(6,'(2A,I4,A,F10.8,A,F10.8,A)') ' NEW SHAKE DID NOT ', & 'CONVERGE AFTER ', i, ' ITERATIONS. ERRX=', & dnrm2(size(sigma), sigma, 1), & ' TOLX=', threshold, '. STOP!' if (.not. cntl%pbicgstab) then write(6,'(A)') " IF PCG IS NOT CONVERGING:" write(6,'(A)') " 1) TRY INCREASING SHAKE_CG_ITER IN THE &CPMD SECTION. EXAMPLE:" write(6,'(A)') " SHAKE_CG_ITER" write(6,'(A)') " 200" write(6,'(A)') " 2) TRY USING THE ALTERNATIVE PBICGSTAB SOLVER. & IN THE &CPMD SECTION, SELECT:" write(6,'(A)') " NEW CONSTRAINTS PBICGSTAB " else write(6,'(A)') " IF PBICGSTAB IS NOT CONVERGING:" write(6,'(A)') " TRY INCREASING SHAKE_CG_ITER IN THE &CPMD SECTION. EXAMPLE:" write(6,'(A)') " SHAKE_CG_ITER" write(6,'(A)') " 200" end if call stopgm('NEW_SHAKE',' ',__LINE__,__FILE__) end if exit end if end do tau0(:,:,:) = taup(:,:,:) call tihalt(procedureN,isub) end subroutine new_shake !> Perform velocity update (within RATTLE algorithm) subroutine new_rattle(matrix, xlagr, tau0, velp, dt, dtm) !> Constraint matrix class(constraint_matrix), intent(inout) :: matrix !> Lagrange multiplier array real(real_8), dimension(:) :: xlagr !> CPMD coordinate matrix real(real_8), dimension(:,:,:), intent(in) :: tau0 !> CPMD velocity matrix real(real_8), dimension(:,:,:), intent(inout) :: velp !> Ionic timestep real(real_8), intent(in) :: dt real(real_8), dimension(:), intent(in) :: dtm real(real_8), dimension(:,:,:), allocatable :: tscr real(real_8), dimension(:), allocatable :: sigma real(real_8), dimension(:, :), allocatable :: grad_sigma real(real_8), dimension(:), allocatable :: dv integer :: i, j, at, atid, spid, cat, pos character(*), PARAMETER :: procedureN = 'new_rattle' integer :: isub call tiset(procedureN,isub) allocate(sigma(size(matrix%nz_row)-1)) allocate(grad_sigma(6, size(matrix%nz_row)-1)) allocate(dv(size(matrix%nz_row)-1)) allocate(tscr(size(velp, 1), size(velp, 2), size(velp, 3))) sigma = 0.0_real_8 grad_sigma = 0.0_real_8 dv = 0.0_real_8 tscr = 0.0_real_8 call dum2(tau0, tscr) call update_constraints(matrix, grad_sigma, sigma, tscr) call fill_matrix(matrix, grad_sigma, dt, dtm) matrix%vals = matrix%vals / dt !\$OMP PARALLEL DO PRIVATE(i, at, cat, spid, atid, pos) do i = 1, size(matrix%nz_row) - 1 do cat = 1, size(matrix%atom_id, 1) at = matrix%atom_id(cat, i) spid = matrix%tau_map(1, cat, i) atid = matrix%tau_map(2, cat, i) pos = (cat - 1) * 3 + 1 dv(i) = dv(i) + dot_product(grad_sigma(pos:pos+2, i), velp(:, atid, spid)) end do end do !\$OMP END PARALLEL DO call pcg_solve(matrix, xlagr, dv, 0.0000001_real_8) dv = 0.0_real_8 !\$OMP PARALLEL DO PRIVATE(i, at, cat, spid, atid, pos) do i = 1, size(matrix%nz_row) - 1 do cat = 1, size(matrix%atom_id, 1) at = matrix%atom_id(cat, i) spid = matrix%tau_map(1, cat, i) atid = matrix%tau_map(2, cat, i) pos = (cat - 1) * 3 + 1 velp(:, atid, spid) = velp(:, atid, spid) - dtm(spid) * xlagr(i) * grad_sigma(pos : pos + 2, i) end do end do !\$OMP END PARALLEL DO call tihalt(procedureN,isub) end subroutine new_rattle !> Use the preconditioned conjugate gradient solver subroutine pcg_solve(matrix, x, b, tol) !> Matrix to solve class(constraint_matrix), intent(in) :: matrix !> Initial guss real(real_8), dimension(:), intent(inout) :: x !> right-hand side real(real_8), dimension(:), intent(in) :: b !> desired tolerance real(real_8) :: tol integer :: i, j real(real_8) :: denominator real(real_8), dimension(:), allocatable :: r real(real_8), dimension(:), allocatable :: precond real(real_8), dimension(:), allocatable :: d, s real(real_8), dimension(:), allocatable :: temp real(real_8) :: delta_start, delta_new, delta_old real(real_8) :: alpha, beta real(real_8), external :: ddot, dnrm2 character(*), PARAMETER :: procedureN = 'pcg_solve' integer :: isub call tiset(procedureN,isub) allocate(r(size(matrix%nz_row)-1)) allocate(precond(size(matrix%nz_row)-1)) allocate(d(size(matrix%nz_row)-1)) allocate(s(size(matrix%nz_row)-1)) allocate(temp(size(matrix%nz_row)-1)) r = 0.0_real_8 precond = 0.0_real_8 d = 0.0_real_8 s = 0.0_real_8 temp = 0.0_real_8 ! Initialize preconditioner which is inverse diagonal !\$OMP PARALLEL DO PRIVATE(i, j) do i = 1, size(matrix%nz_row) - 1 do j = matrix%nz_row(i - 1) + 1, matrix%nz_row(i) if (matrix%r_index(j) == i) then precond(i) = 1.0_real_8 / matrix%vals(j) end if end do end do !\$OMP END PARALLEL DO ! Fill residual and intialize CG variables r = b - matrix * x d = precond * r delta_new = ddot(size(r), r, 1, d, 1) delta_start = delta_new do i = 1, cnti%shake_cg_iter ! Check if we have solved constraints if (dnrm2(size(r), r, 1) < tol) then exit end if temp = matrix * d denominator = ddot(size(d), d, 1, temp, 1) alpha = delta_new / denominator call daxpy(size(d), alpha, d(1), 1, x(1), 1) call daxpy(size(r), -alpha, temp(1), 1, r(1), 1) s = precond * r delta_old = delta_new delta_new = ddot(size(r), r, 1, s, 1) beta = delta_new / delta_old call dscal(size(d), beta, d(1), 1) call daxpy(size(d), 1.0_real_8, s(1), 1, d(1), 1) if (i == cnti%shake_cg_iter) then if (paral%io_parent) then write(6,'(A52,I6,A16)') "WARNING! PCG DID NOT CONVERGE AFTER SHAKE_CG_ITER=", cnti%shake_cg_iter, "STEPS!" end if exit end if end do deallocate(r) deallocate(precond) deallocate(d) deallocate(s) deallocate(temp) call tihalt(procedureN,isub) end subroutine pcg_solve !> Use the preconditioned biconjugate gradient stabilized solver subroutine pbicgstab_solve(matrix, x, b, tol) !> Matrix to solve class(constraint_matrix), intent(in) :: matrix !> Initial guss real(real_8), dimension(:), intent(inout) :: x !> right-hand side real(real_8), dimension(:), intent(in) :: b !> desired tolerance real(real_8) :: tol integer :: i, j real(real_8), dimension(:), allocatable :: precond real(real_8), dimension(:), allocatable :: r, r0, p, y, v, h, s, z, t, temp real(real_8) :: rho_n, rho, alpha, omega, beta real(real_8), external :: ddot, dnrm2 character(*), PARAMETER :: procedureN = 'PBICGSTAB_solve' integer :: isub call tiset(procedureN,isub) ! Allocate and initialize allocate(precond(size(matrix%nz_row)-1)) allocate(r(size(matrix%nz_row)-1)) allocate(r0(size(matrix%nz_row)-1)) allocate(p(size(matrix%nz_row)-1)) allocate(y(size(matrix%nz_row)-1)) allocate(v(size(matrix%nz_row)-1)) allocate(h(size(matrix%nz_row)-1)) allocate(s(size(matrix%nz_row)-1)) allocate(z(size(matrix%nz_row)-1)) allocate(t(size(matrix%nz_row)-1)) allocate(temp(size(matrix%nz_row)-1)) precond = 0.0_real_8 temp = 0.0_real_8 r = 0.0_real_8 r0 = 0.0_real_8 p = 0.0_real_8 y = 0.0_real_8 v = 0.0_real_8 h = 0.0_real_8 s = 0.0_real_8 z = 0.0_real_8 t = 0.0_real_8 beta = 0.0_real_8 rho = 1.0_real_8 rho_n = 1.0_real_8 alpha = 1.0_real_8 omega = 1.0_real_8 ! Initialize preconditioner which is inverse diagonal !\$OMP PARALLEL DO PRIVATE(i, j) do i = 1, size(matrix%nz_row) - 1 do j = matrix%nz_row(i - 1) + 1, matrix%nz_row(i) if (matrix%r_index(j) == i) then precond(i) = 1.0_real_8 / matrix%vals(j) end if end do end do !\$OMP END PARALLEL DO ! Fill residual and intialize CG variables r = b - matrix * x r0 = r ! Check if we have already solved constraints if (dnrm2(size(r), r, 1) < tol) then deallocate(precond) deallocate(temp) deallocate(r) deallocate(r0) deallocate(v) deallocate(y) deallocate(h) deallocate(t) deallocate(s) deallocate(z) deallocate(p) call tihalt(procedureN,isub) return end if ! Main loop do i = 1, cnti%shake_cg_iter rho = rho_n rho_n = ddot(size(r0), r0, 1, r, 1) beta = ( rho_n / rho ) * ( alpha / omega ) p = r + beta * ( p - omega * v ) y = precond * p v = matrix * y alpha = rho_n / ddot(size(r0), r0, 1, v, 1) h = x + alpha * y ! Check convergence temp = b - matrix * h if (dnrm2(size(temp), temp, 1) < tol) then x = h exit end if s = r - alpha * v z = precond * s t = matrix * z omega = ddot(size(t), precond * t, 1, s, 1) / ddot(size(t), precond * t, 1, t, 1) x = h + omega * z ! Check convergence temp = b - matrix * x if (dnrm2(size(temp), temp, 1) < tol) then exit end if r = s - omega * t if (i == cnti%shake_cg_iter) then if (paral%io_parent) then write(6,'(A52,I6,A16)') "WARNING! PBICGSTAB DID NOT CONVERGE AFTER SHAKE_CG_ITER=", cnti%shake_cg_iter, " STEPS!" end if exit end if end do deallocate(precond) deallocate(r) deallocate(r0) deallocate(v) deallocate(y) deallocate(h) deallocate(t) deallocate(s) deallocate(z) deallocate(p) deallocate(temp) call tihalt(procedureN,isub) end subroutine pbicgstab_solve end module constraint_utils . EOF cat << EOF > patch.07 diff ./src/control_def_utils.mod.F90 ./src/control_def_utils.mod.F90 667a cntl%mimic = .FALSE. cntl%new_constraints = .FALSE. cntl%pbicgstab = .FALSE. cnti%shake_maxstep = 5000 cnti%shake_cg_iter = 50 cntl%anneal_dual = .FALSE. cntr%anneal_factors = 1.0_real_8 . EOF cat << EOF > patch.08 diff ./src/control_pri_utils.mod.F90 ./src/control_pri_utils.mod.F90 516,517c IF (cntl%anneal_dual) THEN WRITE(6,'(A,/,A,T54,F12.6,/,A,T54,F12.6)') & ' SIMULATED ANNEALING OF QM/MM WITH DUAL FACTORS', & ' SCALING FACTOR FOR QM SUBSYSTEM: ', cntr%anneal_factors(1), & ' SCALING FACTOR FOR MM SUBSYSTEM: ', cntr%anneal_factors(2) ELSE WRITE(6,'(A,T46,A8,F12.6)')& ' SIMULATED ANNEALING OF IONS WITH ','ANNERI =',cntr%anneri END IF . EOF cat << EOF > patch.09 diff ./src/control_utils.mod.F90 ./src/control_utils.mod.F90 4036a ELSEIF ( keyword_contains(line,'MIMIC') ) THEN #ifdef __MIMIC cntl%mimic = .TRUE. #else something_went_wrong = .true. error_message = 'MIMIC IS REQUESTED BUT CPMD IS NOT COMPILED WITH IT!' #endif ELSEIF ( keyword_contains(line,'NEW', and='CONSTRAINTS')) THEN cntl%new_constraints = .TRUE. IF ( keyword_contains(line,'PBICGSTAB') ) cntl%pbicgstab = .TRUE. !Use pbicgstab (default is pcg) ELSEIF ( keyword_contains(line,'SHAKE_MAXSTEP',alias='SHAKE_MAXSTEPS')) THEN ! Maximum number of shake iterations (cnti%shake_maxstep) previous_line = line READ(iunit,'(A)',iostat=ierr) line CALL readsi(line ,1,last,cnti%shake_maxstep,erread) IF (erread) THEN error_message = "ERROR WHILE READING VALUE" something_went_wrong = .true. go_on_reading = .false. ENDIF ELSEIF ( keyword_contains(line,'SHAKE_CG_ITER')) THEN ! Maximum number of pcg steps in each shake iteration (cnti%shake_cg_iter) previous_line = line READ(iunit,'(A)',iostat=ierr) line CALL readsi(line ,1,last,cnti%shake_cg_iter,erread) IF (erread) THEN error_message = "ERROR WHILE READING VALUE" something_went_wrong = .true. go_on_reading = .false. ENDIF . 3162,3163c cntl%annei=.TRUE. IF ( keyword_contains(previous_line,'DUAL') ) THEN cntl%anneal_dual = .TRUE. cntr%anneal_factors = 1.0_real_8 first = 1 DO i = 1, 2 CALL readsr(line,first,last,cntr%anneal_factors(i),erread) first = last IF (erread) THEN error_message = "ERROR WHILE READING VALUE" something_went_wrong = .true. go_on_reading = .false. ENDIF END DO ELSE CALL readsr(line,1,last,cntr%anneri,erread) END IF . EOF cat << EOF > patch.10 diff ./src/cotr.mod.F90 ./src/cotr.mod.F90 72a INTEGER, ALLOCATABLE, SAVE :: const_nconn(:) INTEGER, ALLOCATABLE, SAVE :: const_conn(:,:) . EOF cat << EOF > patch.11 diff ./src/cpmd.F90 ./src/cpmd.F90 413c IF (cntl%mimic) THEN CALL mimic_ifc_destroy END IF . 163a IF (cntl%mimic) THEN CALL mimic_read_input CALL mimic_ifc_handshake CALL mimic_ifc_request_sizes endif . 99a USE mimic_wrapper, ONLY: mimic_read_input, mimic_ifc_handshake,& mimic_ifc_request_sizes, mimic_ifc_destroy . EOF cat << EOF > patch.12 diff ./src/detdof_utils.mod.F90 ./src/detdof_utils.mod.F90 332c IF ((cotc0%mcnstr.GT.0.OR.cotr007%mrestr.GT.0).AND..NOT.cntl%new_constraints) THEN . 209,234c IF (.NOT.cntl%new_constraints) THEN ALLOCATE(askel(cotc0%nodim,cotc0%mcnstr,12),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) ALLOCATE(mm_askel(cotc0%mcnstr,12),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) ALLOCATE(fc(cotc0%nodim),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) ALLOCATE(fv(cotc0%nodim),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) ALLOCATE(csigm(cotc0%nodim),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) CALL zeroing(askel)!,12*cotc0%nodim*cotc0%mcnstr) CALL zeroing(mm_askel)!,12*cotc0%mcnstr) DO i=1,cotc0%mcnstr ityp=ntcnst(1,i) IF (ityp.EQ.1) THEN csigm(i)=dsigma ELSEIF (ityp.EQ.2) THEN csigm(i)=bsigma ELSEIF (ityp.EQ.3) THEN csigm(i)=tsigma ELSEIF (ityp.EQ.4 .OR. ityp.EQ.10 .OR. ityp.EQ.7) THEN csigm(i)=dsigma ELSEIF(ityp.EQ.5 .OR. ityp.EQ.6 .OR. ityp.EQ.8& .OR. ityp.EQ.9) THEN csigm(i)=tsigma ENDIF IF (ityp .NE. 6 .AND. ityp .NE. 8 .AND.& ityp .NE. 9 .AND. ityp .NE. 10) THEN ia=ntcnst(2,i) ib=ntcnst(3,i) ic=ntcnst(4,i) id=ntcnst(5,i) CALL builda(ia,i,1) CALL builda(ib,i,2) CALL builda(ic,i,3) CALL builda(id,i,4) ENDIF ENDDO ENDIF . 203,206d 180,196d 170c IF ((cotc0%mcnstr.GT.0 .OR. cotr007%mrestr.GT.0).AND..NOT.cntl%new_constraints) THEN . 94c IF (isos1%tisos.AND.(.NOT.lqmmm%qmmm .AND..NOT.cntl%mimic)) THEN . EOF cat << EOF > patch.13 diff ./src/eextern_utils.mod.F90 ./src/eextern_utils.mod.F90 64c IF (.NOT.lqmmm%qmmm.AND..NOT.cntl%mimic)THEN . 60a IF (cntl%mimic.AND.mimic_control%update_potential) THEN mimic_control%update_potential = .FALSE. CALL mimic_ifc_compute_energy() END IF IF (cntl%mimic) THEN IF (tfor) THEN IF (paral%io_parent) THEN CALL mimic_ifc_collect_forces() IF (.NOT.mimic_control%tot_scf_energy) THEN CALL mimic_ifc_collect_energies() END IF END IF CALL mimic_ifc_compute_forces() END IF END IF . 58c IF (tfor.OR.cntl%texadd.OR.cntl%mimic) THEN IF (cntl%mimic) THEN IF (mimic_control%update_potential) THEN CALL mimic_ifc_compute_potential() IF (.NOT.mimic_control%tot_scf_energy) THEN mimic_energy%mm_energy = 0.0_real_8 END IF END IF END IF . 13a USE mimic_wrapper, ONLY: mimic_control,& mimic_energy,& mimic_ifc_collect_energies,& mimic_ifc_collect_forces,& mimic_ifc_compute_energy,& mimic_ifc_compute_forces,& mimic_ifc_compute_potential . EOF cat << EOF > patch.14 diff ./src/finalp_utils.mod.F90 ./src/finalp_utils.mod.F90 77c IF (.NOT.cntl%new_constraints) CALL cnstpr . EOF cat << EOF > patch.15 diff ./src/forcedr_driver.mod.F90 ./src/forcedr_driver.mod.F90 143c CALL mp_bcast(fion,size(fion),parai%io_source,parai%cp_grp) . EOF cat << EOF > patch.16 diff ./src/forces_diag_utils.mod.F90 ./src/forces_diag_utils.mod.F90 368c CALL mp_bcast(fion,size(fion),parai%source,parai%cp_grp) . EOF cat << EOF > patch.17 diff ./src/forces_driver.mod.F90 ./src/forces_driver.mod.F90 425c CALL mp_sum(fion,size(fion),parai%allgrp) . 236a IF (cntl%mimic) THEN mimic_energy%qm_energy = ener_com%etot ener_com%etot = ener_com%etot & + ener_com%eext & + mimic_energy%qmmm_energy & + mimic_energy%mm_energy END IF . 26a USE mimic_wrapper, ONLY: mimic_energy . EOF cat << EOF > patch.18 diff ./src/geofile_utils.mod.F90 ./src/geofile_utils.mod.F90 241a IF (cntl%mimic) THEN CALL mimic_revert_dim() ENDIF . 65a IF (cntl%mimic) THEN CALL mimic_save_dim() CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 15a USE mimic_wrapper, ONLY: mimic_revert_dim,& mimic_save_dim,& mimic_switch_dim . EOF cat << EOF > patch.19 diff ./src/initrun_driver.mod.F90 ./src/initrun_driver.mod.F90 347a IF (cntl%mimic) THEN CALL mimic_revert_dim() ENDIF . 319a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 288a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 202a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 198,199c CALL mp_bcast(tau0,size(tau0),parai%io_source,parai%cp_grp) CALL mp_bcast(velp,size(velp),parai%io_source,parai%cp_grp) . 193a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 190a IF (cntl%mimic) THEN CALL mimic_update_coords(tau0, c0, cm, nstate, ncpw%ngw, inyh) IF (mimic_control%long_range_coupling) THEN CALL mimic_ifc_sort_fragments() END IF CALL mimic_switch_dim(go_qm=.TRUE.) END IF . 178a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 175,176c CALL mp_bcast(tau0,size(tau0),parai%io_source,parai%cp_grp) CALL mp_bcast(velp,size(velp),parai%io_source,parai%cp_grp) . 173a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 172a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 169,170c IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF CALL mp_bcast(tau0,size(tau0),parai%io_source,parai%cp_grp) CALL mp_bcast(velp,size(velp),parai%io_source,parai%cp_grp) IF (cntl%mimic.AND.mimic_control%long_range_coupling) THEN CALL mimic_ifc_sort_fragments() END IF . 164a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 163a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 156a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 155a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 153c IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF CALL mp_bcast(tau0,size(tau0),parai%io_source,parai%cp_grp) IF (cntl%mimic.AND.mimic_control%long_range_coupling) THEN CALL mimic_ifc_sort_fragments() END IF . 148a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 115a IF (cntl%mimic) THEN CALL mimic_save_dim() ENDIF . 25a USE mimic_wrapper, ONLY: mimic_ifc_sort_fragments,& mimic_revert_dim,& mimic_save_dim,& mimic_switch_dim,& mimic_update_coords,& mimic_control . 15a USE cppt, ONLY: inyh . EOF cat << EOF > patch.20 diff ./src/loadpa_utils.mod.F90 ./src/loadpa_utils.mod.F90 432a IF (cntl%mimic) THEN CALL mimic_revert_dim() ENDIF . 117a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 106a IF (cntl%mimic) THEN CALL mimic_save_dim() CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 30c parm, spar, cntl . 18a USE mimic_wrapper, ONLY: mimic_save_dim, & mimic_switch_dim, & mimic_revert_dim . EOF cat << EOF > patch.21 diff ./src/md_driver.mod.F90 ./src/md_driver.mod.F90 1689a IF (cntl%mimic) THEN CALL mimic_revert_dim() ENDIF . 1468a IF (cntl%mimic) THEN CALL mimic_update_coords(taup, c0, cm, nstate, ncpw%ngw, inyh) IF (paral%io_parent) THEN CALL mimic_ifc_send_coordinates() END IF IF (mimic_control%long_range_coupling) THEN CALL mimic_ifc_sort_fragments() END IF END IF . 1404a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 1394a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 1211a IF (cntl%mimic) THEN CALL mimic_subsystem_temperatures(velp) END IF . 1202c IF (paral%parent.AND..NOT.cntl%mimic) CALL geofile(taup,velp,'WRITE') . 1200c IF (cntl%new_constraints) THEN CALL do_rattle(taup, velp) ELSE CALL rattle(taup,velp) END IF . 1189a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 1177a IF (cntl%mimic) THEN CALL mimic_sum_forces(fion) END IF . 1080c ifcalc = 0 . 1044a IF (cntl%mimic) THEN CALL mimic_update_coords(taup, c0, cm, nstate, ncpw%ngw, inyh) mimic_control%update_potential = .TRUE. IF (paral%io_parent) THEN CALL mimic_ifc_send_coordinates() IF (mimic_control%tot_scf_energy) THEN CALL mimic_ifc_collect_energies() END IF END IF IF (mimic_control%long_range_coupling.AND.(mod(iteropt%nfi-1,mimic_control%update_sorting).EQ.0)) THEN CALL mimic_ifc_sort_fragments() END IF CALL mimic_switch_dim(go_qm=.TRUE.) END IF . 1036c IF (cntl%new_constraints) THEN CALL do_shake(tau0, taup, velp) ELSE IF (cotc0%mcnstr.NE.0) CALL cpmdshake(tau0,taup,velp) END IF . 1009a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 956a IF (cntl%mimic) THEN mimic_energy%qm_energy = ener_com%etot ener_com%etot = ener_com%etot & + ener_com%eext & + mimic_energy%qmmm_energy & + mimic_energy%mm_energy END IF . 931a IF (cntl%mimic) THEN CALL mimic_sum_forces(fion) END IF . 874a IF (cntl%mimic) THEN CALL mimic_update_coords(tau0, c0, cm, nstate, ncpw%ngw, inyh) IF (paral%io_parent) THEN CALL mimic_ifc_send_coordinates() IF (mimic_control%tot_scf_energy) THEN CALL mimic_ifc_collect_energies() END IF END IF END IF . 811a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 788c IF (paral%io_parent) THEN IF (cntl%new_constraints) THEN CALL do_rattle(tau0, velp) ELSE CALL rattle(tau0,velp) END IF END IF . 784c IF (paral%io_parent) THEN IF (cntl%new_constraints) THEN CALL do_rattle(tau0, velp) ELSE CALL rattle(tau0,velp) END IF END IF . 770a IF (paral%io_parent.AND.cntl%new_constraints) THEN WRITE(6, '(1x,a)') "USING NEW CONSTRAINTS SOLVER" IF (cntl%pbicgstab) then WRITE(6, '(1x,a)') "WITH PRECONDITIONED BICONJUGATE GRADIENT STABILIZED" ELSE WRITE(6, '(1x,a)') "WITH PRECONDITIONED CONJUGATE GRADIENT" ENDIF WRITE(6, '(1x,a,t60,i6)') "MAX NUMBER OF SHAKE ITERATIONS:", & cnti%shake_maxstep WRITE(6, '(1x,a,t60,i6)') "MAX NUMBER OF PCG STEPS IN EACH SHAKE ITERATION:", & cnti%shake_cg_iter CALL init_constraints(ntcnst, cnsval, ions0%na, ions1%nsp) IF (.NOT.mimic_control%external_constraints) THEN WRITE(6, '(1x,a)') "EXTERNAL CONSTRAINTS DEFINED THROUGH MIMIC WILL BE IGNORED" ENDIF IF (cntl%mimic) THEN CALL mimic_subsystem_dof() END IF END IF . 758a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 673a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 656a IF (cntl%mimic) THEN CALL mimic_ifc_init(rhoe, extf, tau0) mimic_control%update_potential = .TRUE. END IF . 383,388c IF (.NOT.cntl%mimic) THEN ALLOCATE(taup(3,maxsys%nax,maxsys%nsx),stat=ierr) IF(ierr/=0) CALL stopgm(proceduren,'allocation problem',& __LINE__,__FILE__) ALLOCATE(fion(3,maxsys%nax,maxsys%nsx),stat=ierr) IF(ierr/=0) CALL stopgm(proceduren,'allocation problem',& __LINE__,__FILE__) END IF . 290a ELSE IF (cntl%mimic) THEN ALLOCATE(ch(1,1,1),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) ALLOCATE(extf(fpar%kr1*fpar%kr2s*fpar%kr3s),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) . 251a IF (cntl%mimic) THEN CALL mimic_save_dim() CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 171c USE shake_utils, ONLY: cpmdshake,& init_constraints,& do_shake,& do_rattle . 101a USE mimic_wrapper, ONLY: mimic_ifc_collect_energies,& mimic_ifc_collect_forces,& mimic_ifc_sort_fragments,& mimic_ifc_init,& mimic_revert_dim,& mimic_save_dim,& mimic_ifc_send_coordinates,& mimic_sum_forces,& mimic_switch_dim,& mimic_update_coords,& mimic_control,& mimic_energy,& mimic_subsystem_temperatures,& mimic_subsystem_dof . 42c USE cotr, ONLY: cotc0,& cnsval,& ntcnst USE cppt, ONLY: inyh . EOF cat << EOF > patch.22 diff ./src/mdmain_utils.mod.F90 ./src/mdmain_utils.mod.F90 1246a IF (cntl%mimic) THEN CALL mimic_revert_dim() ENDIF . 1150a IF (cntl%mimic) THEN CALL mimic_update_coords(taup, c0, cm, nstate, ncpw%ngw, inyh) IF (paral%io_parent) THEN CALL mimic_ifc_send_coordinates() END IF IF (mimic_control%long_range_coupling) THEN CALL mimic_ifc_sort_fragments() END IF END IF . 1143a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 1139a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 999a IF (cntl%mimic) then CALL mimic_subsystem_temperatures(velp) ENDIF . 926c IF (paral%parent.AND..NOT.cntl%mimic) CALL geofile(taup,velp,'WRITE') . 910c IF (cotc0%mcnstr.NE.0) THEN IF (cntl%new_constraints) THEN CALL do_rattle(taup, velp) ELSE CALL rattle(taup,velp) END IF END IF . 898a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 882a IF (cntl%mimic) THEN CALL mimic_sum_forces(fion) END IF . 739a IF (cntl%mimic) THEN CALL mimic_update_coords(taup, c0, cm, nstate, ncpw%ngw, inyh) mimic_control%update_potential = .TRUE. IF (paral%io_parent) THEN CALL mimic_ifc_send_coordinates() IF (mimic_control%tot_scf_energy) THEN CALL mimic_ifc_collect_energies() END IF END IF IF (mimic_control%long_range_coupling.AND.(mod(iteropt%nfi-1,mimic_control%update_sorting).EQ.0)) THEN CALL mimic_ifc_sort_fragments() END IF CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 729c IF (cntl%new_constraints) THEN CALL do_shake(tau0, taup, velp) ELSE CALL cpmdshake(tau0,taup,velp) ENDIF . 683a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 571a IF (cntl%mimic) THEN CALL mimic_sum_forces(fion) END IF . 528a IF (cntl%mimic) THEN CALL mimic_update_coords(tau0, c0, cm, nstate, ncpw%ngw, inyh) IF (paral%io_parent) THEN CALL mimic_ifc_send_coordinates() IF (mimic_control%tot_scf_energy) THEN CALL mimic_ifc_collect_energies() END IF END IF END IF . 479a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 453c IF (paral%parent) THEN IF (cntl%new_constraints) THEN CALL do_rattle(tau0, velp) ELSE CALL rattle(tau0,velp) END IF END IF . 449c IF (paral%parent) THEN IF (cntl%new_constraints) THEN CALL do_rattle(tau0, velp) ELSE CALL rattle(tau0,velp) END IF END IF . 437a ! INITIALIZE CONSTRAINTS IF (paral%io_parent.AND.cntl%new_constraints) THEN WRITE(6, '(1x,a)') "USING NEW CONSTRAINTS SOLVER" IF (cntl%pbicgstab) then WRITE(6, '(1x,a)') "WITH PRECONDITIONED BICONJUGATE GRADIENT STABILIZED" ELSE WRITE(6, '(1x,a)') "WITH PRECONDITIONED CONJUGATE GRADIENT" ENDIF WRITE(6, '(1x,a,t60,i6)') "MAX NUMBER OF SHAKE ITERATIONS:", & cnti%shake_maxstep WRITE(6, '(1x,a,t60,i6)') "MAX NUMBER OF PCG STEPS IN EACH SHAKE ITERATION:", & cnti%shake_cg_iter CALL init_constraints(ntcnst, cnsval, ions0%na, ions1%nsp) IF (.NOT.mimic_control%external_constraints) THEN WRITE(6, '(1x,a)') "EXTERNAL CONSTRAINTS DEFINED THROUGH MIMIC WILL BE IGNORED" ENDIF IF (cntl%mimic) THEN CALL mimic_subsystem_dof() END IF END IF . 434a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 386a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 346a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 337a IF (cntl%mimic) THEN CALL mimic_ifc_init(rhoe, extf, tau0) mimic_control%update_potential = .TRUE. END IF . 243,250c ! IF (.NOT. cntl%mimic) THEN ALLOCATE(taup(3,maxsys%nax,maxsys%nsx),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) CALL zeroing(taup)!, 3*maxsys%nax*maxsys%nsx) ALLOCATE(fion(3,maxsys%nax,maxsys%nsx),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) ELSE ALLOCATE(extf(fpar%kr1*fpar%kr2s*fpar%kr3s),STAT=ierr) END IF . 218a IF (cntl%mimic) THEN CALL mimic_save_dim() CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 143c USE shake_utils, ONLY: cpmdshake,& init_constraints,& do_shake,& do_rattle . 86a USE mimic_wrapper, ONLY: mimic_ifc_collect_energies,& mimic_ifc_collect_forces,& mimic_ifc_sort_fragments,& mimic_ifc_init,& mimic_revert_dim,& mimic_save_dim,& mimic_ifc_send_coordinates,& mimic_sum_forces,& mimic_switch_dim,& mimic_update_coords,& mimic_control,& mimic_energy,& mimic_subsystem_temperatures,& mimic_subsystem_dof . 43a USE efld, ONLY: extf . 35c USE cotr, ONLY: cnsval,& cotc0,& ntcnst USE cppt, ONLY: inyh . EOF cat << EOF > patch.23 diff ./src/mimic_wrapper.mod.F90 ./src/mimic_wrapper.mod.F90 0a !--------------------------------------------------------------------------------------------------- ! MODULE: mimic_wrapper !> @author Viacheslav Bolnykh - RWTH Aachen/Cyprus Institute !> @author Jogvan Magnus Haugaard Olsen - DTU (jmho@kemi.dtu.dk) ! ! DESCRIPTION: !> Wrapper module for MiMiC main routines !--------------------------------------------------------------------------------------------------- module mimic_wrapper #ifdef __MIMIC use mimic_constants, only: mimic_bohr, mimic_covrad use mimic_main use mimic_tensors, only: initialize_tensors, & terminate_tensors #endif use cell, only: cell_com use control_pri_utils, only: control_pri use cnst, only: factem use coor, only: tau0, taup, fion, velp, lvelini use cotr, only: cnsval, grate, cnpar, cnsval_dest, lskcor, & cotc0, ntcnst, const_nconn, const_conn use efld, only : textfld, extf use error_handling, only: stopgm use gvec, only: gvec_com use inscan_utils, only: inscan use ions, only: ions0, ions1 use isos, only: isos1, isos3 use kinds, only: real_8 use mm_dimmod, only: clsaabox use mm_extrap, only: cold use mp_interface, only: mp_bcast, mp_sum, mp_sync use parac, only: parai, paral use readsr_utils, only: readsr, readsi, input_string_len, keyword_contains use rmas, only: rmass use system, only: maxsys, fpar, parap, cntl, spar, parm, cnti, nkpt, iatpt use timer, only: tihalt, tiset implicit none !> MiMiC control type, maps on the input file section type :: mimic_control_type !> vectors of the whole simulation box real(real_8), dimension(3,3) :: box !> list of covalent radii (read from external file, per species) real(real_8), dimension(:), allocatable :: covalent_radius !> list of factors to be applied to forces contributions of !> different code (one per code) real(real_8), dimension(:), allocatable :: factors !> paths to working folders of each of the client code character(len=250), dimension(:), allocatable :: paths !> total number of species in the simulation integer :: num_species !> maximum nuber of atoms per species integer :: max_atoms !> number of client codes integer :: num_client !> number of quantum species integer :: num_quantum_species !> maximum number of atoms per quantum species integer :: max_quantum_atoms !> total number of atoms integer :: num_atoms !> total number of quantum atoms integer :: num_quantum_atoms !> switch to turn on the computation of potential logical :: update_potential = .false. !> indicator if the dimension is switched to MM system logical :: dim_mm = .false. !> number of saved dimension status integer :: num_save = 0 !> save dimension status logical, dimension(20) :: dim_save = .false. !> use total energy in SCF logical :: tot_scf_energy = .false. !> use long-range coupling scheme logical :: long_range_coupling = .false. !> distance cutoff (in bohr) defining short- and long-range regions real(real_8) :: cutoff_distance = huge(real_8) !> method used to sort fragments integer :: sorting_method = huge(0) !> how often to update the sorted fragments (default every step) integer :: update_sorting = huge(0) !> multipole expansion order integer :: multipole_order = huge(0) !> include constraints from external programs or not logical :: external_constraints = .true. !> degrees of freedom for QM subsystem real(real_8) :: qm_dof !> degrees of freedom for MM subsystem real(real_8) :: mm_dof end type mimic_control_type !> energy holder for QM/MM simulation type :: mimic_energy_type !> QM/MM energy real(real_8) :: qmmm_energy = 0.0_real_8 !> QM energy real(real_8) :: qm_energy = 0.0_real_8 !> MM energy real(real_8) :: mm_energy = 0.0_real_8 end type mimic_energy_type ! IDs of the first and the last SR atom of the CP group integer, dimension(:), allocatable, private :: gr_sr_atom_start, gr_sr_atom_end ! IDs of the first and the last LR atom of the CP group integer, dimension(:), allocatable, private :: gr_lr_atom_start, gr_lr_atom_end ! IDs of the first and the last SR atom of the current process integer, dimension(:), allocatable, private :: pr_sr_atom_start, pr_sr_atom_end ! IDs of the first and the last LR atom of the current process integer, dimension(:), allocatable, private :: pr_lr_atom_start, pr_lr_atom_end !> internal forces array to store QM/MM and MM forces real(real_8), dimension(:,:,:), allocatable :: mimic_forces !> internal forces array to store MM forces real(real_8), dimension(:,:,:), allocatable :: mm_forces !> coordinates array storing all of the coordinates with MM box PBC applied real(real_8), dimension(:,:,:), allocatable :: wrapped_coordinates !> MiMiC options type(mimic_control_type), save :: mimic_control !> storage for the energy type(mimic_energy_type), save :: mimic_energy #ifdef __MIMIC !> temporary storage for the allocation sizes type(sizes_type), save :: sizes !> system information type(system_type), save :: system_data !> array of subsystems associated with each client code type(subsystem_type), dimension(:), allocatable :: subsystems !> quantum subsystem type(quantum_fragment_type) :: quantum_fragment #else character(len=*), parameter :: MIMIC_MISSING_ERROR = "CPMD is compiled without MiMiC support but MiMiC procedures are called!" #endif procedure(cpmd_error_handler), pointer, save :: error_handler => cpmd_error_handler contains !> subroutine to handle MiMiC errors subroutine cpmd_error_handler(err_type, message, source_file, line_num) !> Type of the error integer, intent(in) :: err_type !> Optional message character(len=*), intent(in) :: message !> Source file where the problem occurred character(len=*), intent(in) :: source_file !> Source line at which the error occurred integer, intent(in) :: line_num call stopgm("cpmd_error_handler", & "MiMiC has encountered an error, see MiMiC_error file(s) for details", & line_num, source_file) end subroutine cpmd_error_handler !> initialize overlaps within MiMiC subroutine mimic_ifc_init_overlaps(overlaps) !> overlap maps integer, dimension(:,:), intent(in) :: overlaps character(*), parameter :: procedureN = 'mimic_ifc_init_overlaps' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC call mimic_init_overlaps(overlaps) #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_ifc_init_overlaps !> finalize MiMiC simulation and destroy internal objects subroutine mimic_ifc_destroy character(*), parameter :: procedureN = 'mimic_ifc_destroy' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC call terminate_tensors() if (paral%io_parent) then call mimic_finalize call mimic_destroy end if #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_ifc_destroy #ifdef __MIMIC !> request system information: coordinates, multipoles, etc. subroutine mimic_ifc_request_system_data(sizes, system_data) !> sizes of data structures type(sizes_type), intent(inout) :: sizes !> system information type(system_type), intent(inout) :: system_data character(*), parameter :: procedureN = 'mimic_ifc_request_system_data' integer :: isub call tiset(procedureN, isub) call mimic_request_system_data(sizes, system_data) call tihalt(procedureN, isub) end subroutine mimic_ifc_request_system_data #endif !> collect forces from clients subroutine mimic_ifc_collect_forces() character(*), parameter :: procedureN = 'mimic_ifc_collect_forces' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC call mimic_collect_forces(subsystems) #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_ifc_collect_forces !> collect energies from clients subroutine mimic_ifc_collect_energies() character(*), parameter :: procedureN = 'mimic_ifc_collect_energies' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC call mimic_collect_energies(subsystems, mimic_energy%mm_energy) #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_ifc_collect_energies !> compute QM/MM interaction energy subroutine mimic_ifc_compute_energy real(real_8) :: temp_energy real(real_8), dimension(:,:), allocatable :: tensors, tensor_sums character(*), parameter :: procedureN = 'mimic_ifc_compute_energy' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC if (mimic_control%long_range_coupling) then call mimic_compute_multipoles(quantum_fragment, mimic_control%multipole_order, & parap%NRXPL(parai%mepos,1), parap%NRXPL(parai%mepos,2)) call mp_sum(quantum_fragment%electronic_multipoles%values, & size(quantum_fragment%electronic_multipoles%values), & parai%allgrp) call mimic_compute_tensors(subsystems, & quantum_fragment, & tensors, & tensor_sums, & pr_lr_atom_start, pr_lr_atom_end) end if mimic_energy%qmmm_energy = 0.0_real_8 call mimic_compute_energy(subsystems, & quantum_fragment, & mimic_energy%qmmm_energy, & pr_sr_atom_start, & pr_sr_atom_end, & tensor_sums) CALL mp_sum(mimic_energy%qmmm_energy, parai%cp_grp) #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_ifc_compute_energy !> compute QM/MM forces subroutine mimic_ifc_compute_forces() character(*), parameter :: procedureN = 'mimic_ifc_compute_forces' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC if (mimic_control%long_range_coupling) then call mimic_compute_multipoles(quantum_fragment, mimic_control%multipole_order, & parap%NRXPL(parai%mepos,1), parap%NRXPL(parai%mepos,2)) call mp_sum(quantum_fragment%electronic_multipoles%values, & size(quantum_fragment%electronic_multipoles%values), & parai%allgrp) end if call mimic_compute_forces(subsystems, & quantum_fragment, & parap%NRXPL(parai%mepos,1), parap%NRXPL(parai%mepos,2), & gr_sr_atom_start, gr_sr_atom_end, & gr_lr_atom_start, gr_lr_atom_end, & pr_sr_atom_start, pr_sr_atom_end, & pr_lr_atom_start, pr_lr_atom_end, paral%parent) #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_ifc_compute_forces !> compute external potential due to MM atoms subroutine mimic_ifc_compute_potential() real(real_8), dimension(:,:), allocatable :: tensors, tensor_sums character(*), parameter :: procedureN = 'mimic_ifc_compute_potential' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC if (mimic_control%long_range_coupling) then call mimic_compute_tensors(subsystems, & quantum_fragment, & tensors, & tensor_sums, & pr_lr_atom_start, pr_lr_atom_end) call mp_sum(tensor_sums, size(tensor_sums), parai%allgrp) tensor_sums = -tensor_sums end if call mimic_compute_potential(subsystems, & quantum_fragment, & parap%NRXPL(parai%mepos,1), parap%NRXPL(parai%mepos,2), & gr_sr_atom_start, gr_sr_atom_end, tensor_sums) CALL mp_sum(extf, size(extf), parai%cp_inter_grp) #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_ifc_compute_potential !> sort fragments into short- and long-ranged shells subroutine mimic_ifc_sort_fragments() integer :: i integer :: tot_sr_atoms integer :: tot_lr_atoms integer :: num_sr_atoms integer :: num_lr_atoms integer :: atom_start integer :: atom_end integer :: remainder character(*), parameter :: procedureN = 'mimic_ifc_sort_fragments' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC tot_sr_atoms = 0 tot_lr_atoms = 0 if (paral%io_parent) then write(6,'(8x,a)') 'SORTING MM ATOMS INTO SHORT- AND LONG-RANGE DOMAINS' end if do i = 1, size(subsystems) call mimic_sort_fragments(subsystems(i), quantum_fragment, & mimic_control%cutoff_distance, mimic_control%sorting_method) call mimic_distribute_atoms(subsystems(i)%num_sr_atoms, parai%cp_nogrp, & parai%cp_inter_me, gr_sr_atom_start(i), gr_sr_atom_end(i)) call mimic_distribute_atoms(subsystems(i)%num_lr_atoms, parai%cp_nogrp, & parai%cp_inter_me, gr_lr_atom_start(i), gr_lr_atom_end(i)) call mimic_distribute_atoms(subsystems(i)%num_sr_atoms, parai%cp_nproc, & parai%cp_me, pr_sr_atom_start(i), pr_sr_atom_end(i)) call mimic_distribute_atoms(subsystems(i)%num_lr_atoms, parai%cp_nproc, & parai%cp_me, pr_lr_atom_start(i), pr_lr_atom_end(i)) if (paral%io_parent) then ! write(6,'(8x,a,i2,a)') "SUBSYSTEM ", i, ":" write(6,'(8x,a,t58,i8)') "SHORT-RANGE ATOMS PER PROCESS GROUP:", & gr_sr_atom_end(i) - gr_sr_atom_start(i) + 1 write(6,'(8x,a,t58,i8)') "SHORT-RANGE ATOMS PER PROCESS:", & pr_sr_atom_end(i) - pr_sr_atom_start(i) + 1 write(6,'(8x,a,t58,i8)') "LONG-RANGE ATOMS PER PROCESS GROUP:", & gr_lr_atom_end(i) - gr_lr_atom_start(i) + 1 write(6,'(8x,a,t58,i8)') "LONG-RANGE ATOMS PER PROCESS:", & pr_lr_atom_end(i) - pr_lr_atom_start(i) + 1 end if end do #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_ifc_sort_fragments !> Distribute atoms across the processes of the given communicator subroutine mimic_distribute_atoms(tot_num_atoms, comm_size, proc_id, atom_start, atom_end) !> total number of atoms to distribute integer, intent(in) :: tot_num_atoms !> size of the communicator, i.e., number of processes in the communicator integer, intent(in) :: comm_size !> ID of the current process (starts counting from zero) integer, intent(in) :: proc_id !> index marking the start of the chunk of the atom array to be treated by this process integer, intent(out) :: atom_start !> index marking the end of the chunk of the atom array to be treated by this process integer, intent(out) :: atom_end integer :: remainder integer :: n_atoms integer :: temp_final character(*), parameter :: procedureN = 'mimic_distribute_atoms' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC if (proc_id > comm_size) then CALL stopgm(procedureN, 'process id does not belong to the communicator', __LINE__, __FILE__) endif n_atoms = tot_num_atoms / comm_size remainder = modulo(tot_num_atoms, comm_size) atom_start = 0 if (remainder > 0) then if (proc_id < remainder) then n_atoms = n_atoms + 1 else atom_start = remainder endif end if atom_start = atom_start + proc_id * n_atoms + 1 temp_final = atom_start + n_atoms - 1 if (temp_final > tot_num_atoms) then temp_final = tot_num_atoms endif atom_end = temp_final #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_distribute_atoms !> add internal forces array to global forces subroutine mimic_sum_forces(fion) !> global forces real(real_8), dimension(:,:,:), intent(inout) :: fion real(real_8), dimension(:,:,:), allocatable :: fion_temp character(*), parameter :: procedureN = 'mimic_sum_forces' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC allocate(fion_temp, mold=mimic_forces) CALL mp_sum(mimic_forces, fion_temp, size(mimic_forces), parai%cp_grp) fion = fion + fion_temp mimic_forces = 0.0_real_8 #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_sum_forces !> save dimensions of atomic data structures subroutine mimic_save_dim() character(*), parameter :: procedureN = 'mimic_save_dim' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC mimic_control%num_save = mimic_control%num_save + 1 mimic_control%dim_save(mimic_control%num_save) = mimic_control%dim_mm #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_save_dim !> revert any changes to dimensions of data structures subroutine mimic_revert_dim() character(*), parameter :: procedureN = 'mimic_revert_dim' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC call mimic_switch_dim(go_qm= .not. mimic_control%dim_save(mimic_control%num_save)) mimic_control%num_save = mimic_control%num_save - 1 #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_revert_dim !> switch dimensions of atomic data structures subroutine mimic_switch_dim(go_qm) !> flag indicating switch MM->QM or vice versa logical, intent(in) :: go_qm character(*), parameter :: procedureN = 'mimic_switch_dim' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC if (.not. cntl%mimic) return if (go_qm) then maxsys%nax = mimic_control%max_quantum_atoms maxsys%nsx = mimic_control%num_quantum_species ions1%nsp = mimic_control%num_quantum_species ions1%nat = mimic_control%num_quantum_atoms mimic_control%dim_mm = .false. else maxsys%nax = mimic_control%max_atoms maxsys%nsx = mimic_control%num_species ions1%nsp = mimic_control%num_species ions1%nat = mimic_control%num_atoms mimic_control%dim_mm = .true. end if #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_switch_dim !> read &MIMIC section of the input file subroutine mimic_read_input() integer, parameter :: iunit = 5 integer, parameter :: output_unit = 6 integer, parameter :: max_unknown_lines = 250 character(len=input_string_len) :: line, previous_line, error_message, & unknown(max_unknown_lines) integer :: first, last, ierr, num_unknown_lines logical :: error, something_went_wrong, go_on_reading integer :: i, j integer :: num_overlaps, num_factors integer, dimension(:, :), allocatable :: overlaps character(len=*), parameter :: procedureN = 'mimic_read_input' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC if (paral%io_parent) then ! write(output_unit,'(1x,a)') 'READING &MIMIC SECTION' num_unknown_lines = 0 line = ' ' previous_line = ' ' error_message = ' ' ierr = inscan(iunit, '&MIMIC') if (ierr == 0) then textfld = .true. mimic_control%tot_scf_energy = .false. mimic_control%long_range_coupling = .false. mimic_control%cutoff_distance = huge(real_8) mimic_control%sorting_method = CENTROID mimic_control%update_sorting = 1 mimic_control%multipole_order = 2 go_on_reading = .true. something_went_wrong = .false. do while (go_on_reading) previous_line = line read(iunit, '(a80)', iostat=ierr) line if (ierr /= 0) then something_went_wrong = .true. go_on_reading = .false. else if (keyword_contains(line, '&END')) then go_on_reading = .false. else if (keyword_contains(line, 'PATHS')) then read(iunit, '(a)', iostat=ierr) line if (ierr /= 0) then something_went_wrong = .true. go_on_reading = .false. end if first = 1 call readsi(line, first, last, mimic_control%num_client, error) if (error) then error_message = 'ERROR WHILE READING VALUE' something_went_wrong = .true. go_on_reading = .false. end if allocate(mimic_control%paths(mimic_control%num_client), stat=ierr) if (ierr /= 0) call stopgm(procedureN, 'Error while allocating array, cf. output file', __LINE__, __FILE__) allocate(mimic_control%factors(mimic_control%num_client), stat=ierr) if (ierr /= 0) call stopgm(procedureN, 'Error while allocating array, cf. output file', __LINE__, __FILE__) mimic_control%factors(:) = 1 do i = 1, mimic_control%num_client read(iunit, fmt='(a)', iostat=ierr) mimic_control%paths(i) if (ierr /= 0) then something_went_wrong = .true. go_on_reading = .false. end if end do else if (keyword_contains(line, 'DISABLE', and='CONSTRAINTS')) then mimic_control%external_constraints = .false. else if (keyword_contains(line, 'OVERLAPS')) then read(iunit, '(a)', iostat=ierr) line if (ierr /= 0) then something_went_wrong = .true. go_on_reading = .false. end if first = 1 call readsi(line, first, last, num_overlaps, error) if (error) then error_message = 'ERROR WHILE READING VALUE' something_went_wrong = .true. go_on_reading = .false. end if allocate(overlaps(4, num_overlaps), stat=ierr) if (ierr /= 0) call stopgm(procedureN, 'Error while allocating array, cf. output file', __LINE__, __FILE__) do i = 1, num_overlaps read(iunit, '(a80)', iostat=ierr) line if (ierr /= 0) then something_went_wrong = .true. go_on_reading = .false. end if first = 1 do j = 1, 4 call readsi(line, first, last, overlaps(j, i), error) if (error) then error_message = 'ERROR WHILE READING VALUE' something_went_wrong = .true. go_on_reading = .false. end if first = last + 1 end do end do else if (keyword_contains(line, 'BOX')) then read(iunit, '(a)', iostat=ierr) line if (ierr /= 0) then something_went_wrong = .true. go_on_reading = .false. end if first = 1 call readsr(line, first, last, mimic_control%box(1,1), error) if (error) then error_message = 'ERROR WHILE READING VALUE' something_went_wrong = .true. go_on_reading = .false. end if first = last + 1 call readsr(line, first, last, mimic_control%box(2,2), error) if (error) then error_message = 'ERROR WHILE READING VALUE' something_went_wrong = .true. go_on_reading = .false. end if first = last + 1 call readsr(line, first, last, mimic_control%box(3,3), error) if (error) then error_message = 'ERROR WHILE READING VALUE' something_went_wrong = .true. go_on_reading = .false. end if else if (keyword_contains(line, 'SUBSYSTEM', and='FACTORS')) then read(iunit, '(a)', iostat=ierr) line if (ierr /= 0) then something_went_wrong = .true. go_on_reading = .false. end if call readsi(line, first, last, num_factors, error) if (num_factors /= mimic_control%num_client) then error_message = 'NUMBER OF FACTORS SHOULD BE EQUAL TO THE NUMBER OF CLIENTS' something_went_wrong = .true. go_on_reading = .false. endif do i = 1, num_factors read(iunit, '(a80)', iostat=ierr) line if (ierr /= 0) then something_went_wrong = .true. go_on_reading = .false. end if first = 1 call readsr(line, first, last, mimic_control%factors(i), error) end do else if (keyword_contains(line, 'TOTAL', and='SCF')) then mimic_control%tot_scf_energy = .true. else if (keyword_contains(line, 'FRAGMENT', and='SORTING')) then if (keyword_contains(line, 'CENTROID')) then mimic_control%sorting_method = CENTROID else if (keyword_contains(line, 'CENTER-OF-MASS')) then mimic_control%sorting_method = CENTER_OF_MASS else if (keyword_contains(line, 'CENTER-OF-CHARGE')) then mimic_control%sorting_method = CENTER_OF_CHARGE else if (keyword_contains(line, 'MINIMUM', and='DISTANCE')) then mimic_control%sorting_method = MINIMUM_DISTANCE else if (keyword_contains(line, 'ATOM-WISE')) then mimic_control%sorting_method = ATOM_WISE else mimic_control%sorting_method = CENTROID end if if (keyword_contains(line, 'UPDATE')) then read(iunit, '(a)', iostat=ierr) line if (ierr /= 0) then something_went_wrong = .true. go_on_reading = .false. end if first = 1 call readsi(line, first, last, mimic_control%update_sorting, error) if (error) then error_message = 'ERROR WHILE READING VALUE' something_went_wrong = .true. go_on_reading = .false. end if end if else if (keyword_contains(line, 'LONG-RANGE', and='COUPLING')) then mimic_control%long_range_coupling = .true. else if (keyword_contains(line, 'CUTOFF', and='DISTANCE')) then read(iunit, '(a)', iostat=ierr) line if (ierr /= 0) then something_went_wrong = .true. go_on_reading = .false. end if first = 1 call readsr(line, first, last, mimic_control%cutoff_distance, error) if (error) then error_message = 'ERROR WHILE READING VALUE' something_went_wrong = .true. go_on_reading = .false. end if else if (keyword_contains(line, 'MULTIPOLE', and='ORDER')) then read(iunit, '(a)', iostat=ierr) line if (ierr /= 0) then something_went_wrong = .true. go_on_reading = .false. end if first = 1 call readsi(line, first, last, mimic_control%multipole_order, error) if (error) then error_message = 'ERROR WHILE READING VALUE' something_went_wrong = .true. go_on_reading = .false. end if else ! Unknown keyword if (' ' /= line) then if (num_unknown_lines < max_unknown_lines) then num_unknown_lines = num_unknown_lines + 1 unknown(num_unknown_lines) = line else do i = 1, max_unknown_lines - 1 unknown(i) = unknown(i+1) enddo unknown(num_unknown_lines) = line endif endif end if end do else something_went_wrong = .true. error_message = 'MISSING &MIMIC SECTION - SECTION MANDATORY FOR MIMIC' end if if (something_went_wrong) THEN write(output_unit,'(/,1x,64("!"))') write(output_unit,'(1x, a, 1x, a)') 'ERROR:', 'PROBLEM WHILE READING &MIMIC SECTION:' write(output_unit,'(8x,a)') trim(adjustl(error_message)) if (line /= ' ' .or. previous_line /= ' ') THEN write(output_unit,'(8x,a)') 'THE LAST TWO LINES READ WITHIN THE SECTION WERE:' write(output_unit,'(/,1x,a)') trim(adjustl(previous_line)) write(output_unit,'(1x,a)') trim(Adjustl(line)) endif write(output_unit,'(1x,64("!"))') call stopgm(procedureN, 'Error while reading &MIMIC section, cf. output file', __LINE__, __FILE__) endif if (num_unknown_lines > 0) then write(6,'(/,1x,64("="))') write(6,'(1x,a,14x,a,12x,a)') '= ', 'UNKNOWN KEYWORDS IN &MIMIC SECTION', ' =' do i = 1, num_unknown_lines write(6,'(1x,a2,a60,a2)') '= ', unknown(i), ' =' end do write(6,'(1x,64("="),/)') end if ! write(output_unit,'(1x,a)') 'DONE READING &MIMIC SECTION' end if if (.not. mimic_control%long_range_coupling) then mimic_control%cutoff_distance = huge(real_8) end if if (cntl%mimic) then call mimic_init_error_handler(parai%me, error_handler) call mp_bcast(textfld, parai%io_source, parai%cp_grp) call mp_bcast(rmass%pma0, size(rmass%pma0), parai%io_source, parai%cp_grp) call mp_bcast(ions0%iatyp, size(ions0%iatyp), parai%io_source, parai%cp_grp) call mp_bcast(mimic_control%num_client, parai%io_source, parai%cp_grp) if (.not. allocated(mimic_control%factors)) allocate(mimic_control%factors(mimic_control%num_client)) call mp_bcast(mimic_control%factors, size(mimic_control%factors), parai%io_source, parai%cp_grp) call mp_bcast(mimic_control%box, size(mimic_control%box), parai%io_source, parai%cp_grp) call mp_bcast(mimic_control%long_range_coupling, parai%io_source, parai%cp_grp) call mp_bcast(mimic_control%sorting_method, parai%io_source, parai%cp_grp) call mp_bcast(mimic_control%update_sorting, parai%io_source, parai%cp_grp) call mp_bcast(mimic_control%cutoff_distance, parai%io_source, parai%cp_grp) call mp_bcast(mimic_control%multipole_order, parai%io_source, parai%cp_grp) mimic_control%max_quantum_atoms = maxsys%nax mimic_control%num_quantum_species = maxsys%nsx call mp_bcast(num_overlaps, parai%io_source, parai%cp_grp) if (.not. allocated(overlaps)) allocate(overlaps(4, num_overlaps)) call mp_bcast(overlaps, size(overlaps), parai%io_source, parai%cp_grp) call mimic_set_num_clients(mimic_control%num_client) allocate(gr_sr_atom_start(mimic_control%num_client)) allocate(gr_sr_atom_end(mimic_control%num_client)) allocate(gr_lr_atom_start(mimic_control%num_client)) allocate(gr_lr_atom_end(mimic_control%num_client)) allocate(pr_sr_atom_start(mimic_control%num_client)) allocate(pr_sr_atom_end(mimic_control%num_client)) allocate(pr_lr_atom_start(mimic_control%num_client)) allocate(pr_lr_atom_end(mimic_control%num_client)) allocate(subsystems(mimic_control%num_client)) call mimic_ifc_init_overlaps(OVERLAPS) end if #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_read_input !> perform handshake operation with clients subroutine mimic_ifc_handshake character(len=*), parameter :: procedureN = 'mimic_ifc_handshake' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC if (paral%io_parent) then call mimic_handshake(mimic_control%paths) end if #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_ifc_handshake !> request allocation sizes for internal data structures subroutine mimic_ifc_request_sizes integer :: max_atom_pfrag integer, dimension(:), allocatable :: atoms_pcode integer, dimension(2) :: sp_map_shape integer :: i, j integer :: cp_inter_grp integer :: ierr character(len=*), parameter :: procedureN = 'mimic_ifc_request_sizes' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC allocate(atoms_pcode(mimic_control%num_client), stat=ierr) if (ierr /= 0) call stopgm(procedureN, 'Error while allocating array, cf. output file', __LINE__, __FILE__) if (paral%io_parent) then call mimic_request_sizes(sizes, system_data) max_atom_pfrag = size(sizes%atoms_pfragment, dim=1) else allocate(sizes%atoms_pcode(mimic_control%num_client)) allocate(sizes%multipoles_order(mimic_control%num_client)) allocate(sizes%multipoles_patom(mimic_control%num_client)) allocate(sizes%frag_num(mimic_control%num_client)) allocate(sizes%nbonds_pcode(mimic_control%num_client)) allocate(sizes%nangles_pcode(mimic_control%num_client)) allocate(sizes%types_length_pcode(mimic_control%num_client)) sizes%nbonds_pcode(:) = 0 sizes%nangles_pcode(:) = 0 end if call mp_sync(parai%cp_grp) call mp_bcast(sizes%atoms_pcode, & mimic_control%num_client, & parai%io_source, & parai%cp_grp) call mp_bcast(sizes%multipoles_order, & mimic_control%num_client, & parai%io_source, & parai%cp_grp) call mp_bcast(sizes%multipoles_patom, & mimic_control%num_client, & parai%io_source, & parai%cp_grp) call mp_bcast(sizes%frag_num, & mimic_control%num_client, & parai%io_source, & parai%cp_grp) call mp_bcast(sizes%types_length_pcode, & mimic_control%num_client, & parai%io_source, & parai%cp_grp) call mp_bcast(max_atom_pfrag, & parai%io_source, & parai%cp_grp) if (.not. allocated(sizes%atoms_pfragment)) & allocate(sizes%atoms_pfragment(max_atom_pfrag, mimic_control%num_client)) call mp_bcast(sizes%atoms_pfragment, size(sizes%atoms_pfragment), parai%io_source, parai%cp_grp) call mp_bcast(sizes%num_atoms, parai%io_source, parai%cp_grp) call mp_bcast(sizes%num_species, parai%io_source, parai%cp_grp) if (.not. allocated(sizes%atoms_pspecies)) allocate(sizes%atoms_pspecies(sizes%num_species)) call mp_bcast(sizes%atoms_pspecies, sizes%num_species, parai%io_source, parai%cp_grp) call mp_bcast(ions0%na, size(ions0%na), parai%io_source, parai%cp_grp) call mp_bcast(ions1%nat, parai%io_source, parai%cp_grp) call mp_bcast(ions1%nsp, parai%io_source, parai%cp_grp) ions1%nsp = maxsys%nsx mimic_control%num_species = maxsys%nsx + sizes%num_species mimic_control%max_atoms = max(maxsys%nax, sizes%num_atoms) mimic_control%num_quantum_atoms = sum(ions0%na) mimic_control%num_atoms = mimic_control%num_quantum_atoms + sum(sizes%atoms_pspecies) if (allocated(tau0)) deallocate(tau0, stat=ierr) if (allocated(velp)) deallocate(velp, stat=ierr) if (allocated(lvelini)) deallocate(lvelini, stat=ierr) if (allocated(taup)) deallocate(taup, stat=ierr) if (allocated(fion)) deallocate(fion, stat=ierr) allocate(tau0(3, mimic_control%max_atoms, mimic_control%num_species), stat=ierr) allocate(velp(3, mimic_control%max_atoms, mimic_control%num_species), stat=ierr) allocate(lvelini(0:mimic_control%max_atoms + 1, mimic_control%num_species), stat=ierr) allocate(taup(3, mimic_control%max_atoms, mimic_control%num_species), stat=ierr) allocate(fion(3, mimic_control%max_atoms, mimic_control%num_species), stat=ierr) tau0 = 0.0_real_8 velp = 0.0_real_8 lvelini = .false. fion = 0.0_real_8 taup = 0.0_real_8 if ( paral%io_parent ) then call mimic_ifc_request_system_data(sizes, system_data) sp_map_shape = shape(system_data%species_map) else allocate(system_data%species(maxval(sizes%atoms_pcode), & mimic_control%num_client)) allocate(system_data%masses(sizes%num_species)) allocate(system_data%elements(sizes%num_species)) allocate(system_data%multipole_values(maxval(sizes%multipoles_patom), & maxval(sizes%atoms_pcode), & mimic_control%num_client)) allocate(system_data%atom_fragment_ids(maxval(sizes%atoms_pfragment), & maxval(sizes%frag_num), & mimic_control%num_client)) end if call mp_bcast(system_data%species, & size(system_data%species), & parai%io_source, & parai%cp_grp) call mp_bcast(system_data%masses, & size(system_data%masses), & parai%io_source, & parai%cp_grp) call mp_bcast(system_data%elements, & size(system_data%elements), & parai%io_source, & parai%cp_grp) call mp_bcast(system_data%multipole_values, & size(system_data%multipole_values), & parai%io_source, & parai%cp_grp) call mp_bcast(sp_map_shape, & size(sp_map_shape), & parai%io_source, & parai%cp_grp) if (.not. allocated(system_data%species_map)) & allocate(system_data%species_map(sp_map_shape(1), sp_map_shape(2))) call mp_bcast(system_data%species_map, & size(system_data%species_map), & parai%io_source, & parai%cp_grp) do i = 1, mimic_control%num_client call mp_bcast(system_data%atom_fragment_ids(:,:,i), & size(system_data%atom_fragment_ids(:,:,i)), & parai%io_source, parai%cp_grp) end do allocate(mimic_control%covalent_radius(maxsys%nsx + sizes%num_species)) do i = 1, sizes%num_species mimic_control%covalent_radius(maxsys%nsx + i) = mimic_covrad(system_data%elements(i)) * mimic_bohr end do rmass%pma0(maxsys%nsx + 1:maxsys%nsx + sizes%num_species) = system_data%masses ions0%iatyp(maxsys%nsx + 1:maxsys%nsx + sizes%num_species) = system_data%elements ions0%na(maxsys%nsx + 1 : mimic_control%num_species) = sizes%atoms_pspecies #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_ifc_request_sizes !> initialize constraints information subroutine mimic_init_constraints integer, dimension(:, :), allocatable :: temp_cnst real(real_8), dimension(:), allocatable :: temp_cnst_val real(real_8), dimension(:), allocatable :: temp_grate real(real_8), dimension(:), allocatable :: temp_cnsval_dest real(real_8), dimension(:, :), allocatable :: temp_cnpar integer :: tot_const_num integer :: n_code, n_bond, n_angle integer :: offset character(len=*), parameter :: procedureN = 'mimic_ifc_request_sizes' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC tot_const_num = cotc0%mcnstr if (mimic_control%external_constraints) then do n_code = 1, mimic_control%num_client tot_const_num = tot_const_num + & sizes%nbonds_pcode(n_code) + & sizes%nangles_pcode(n_code) end do endif if (tot_const_num == 0) then return end if if (allocated(ntcnst)) then call move_alloc(ntcnst, temp_cnst) call move_alloc(cnsval, temp_cnst_val) call move_alloc(cnpar, temp_cnpar) call move_alloc(cnsval_dest, temp_cnsval_dest) call move_alloc(grate, temp_grate) end if allocate(ntcnst(6, tot_const_num)) allocate(cnsval(tot_const_num)) allocate(cnpar(2, tot_const_num)) allocate(cnsval_dest(tot_const_num)) allocate(grate(tot_const_num)) ntcnst(:,:) = 0.0_real_8 grate(:) = 0.0_real_8 cnpar(:, :) = 0.0_real_8 cnsval_dest(:) = 0.0_real_8 cnsval(:) = 0.0_real_8 if (allocated(temp_cnst))then ntcnst(1:6, 1:cotc0%mcnstr) = temp_cnst(1:6, 1:cotc0%mcnstr) cnsval(1:cotc0%mcnstr) = temp_cnst_val grate(1:cotc0%mcnstr) = temp_grate cnpar(1:2, 1:cotc0%mcnstr) = temp_cnpar cnsval_dest(1:cotc0%mcnstr) = temp_cnsval_dest end if offset = cotc0%mcnstr + 1 if (mimic_control%external_constraints) then do n_code = 1, mimic_control%num_client do n_bond = 1, sizes%nbonds_pcode(n_code) ntcnst(1, offset) = 1 ntcnst(2, offset) = system_data%bonds(n_bond, n_code)%atom_i ntcnst(3, offset) = system_data%bonds(n_bond, n_code)%atom_j cnsval(offset) = system_data%bonds(n_bond, n_code)%length offset = offset + 1 end do do n_angle = 1, sizes%nangles_pcode(n_code) ntcnst(1, offset) = 2 ntcnst(2, offset) = system_data%angles(n_angle, n_code)%atom_i ntcnst(3, offset) = system_data%angles(n_angle, n_code)%atom_j ntcnst(4, offset) = system_data%angles(n_angle, n_code)%atom_k cnsval(offset) = system_data%angles(n_angle, n_code)%angle offset = offset + 1 end do end do endif cotc0%mcnstr = tot_const_num call mimic_count_trcnst(const_nconn, const_conn, system_data%bonds, sizes%nbonds_pcode) #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_init_constraints !> intialize MiMiC's data structures subroutine mimic_ifc_init(rhoe, extf, tau) !> electronic density array real(real_8), dimension(fpar%kr1, fpar%kr2s, fpar%kr3s), intent(in) :: rhoe !> external potential acting on the electronic grid real(real_8), dimension(fpar%kr1, fpar%kr2s, fpar%kr3s), intent(in) :: extf !> CPMD coordinates array real(real_8), dimension(:,:,:), intent(inout) :: tau real(real_8), dimension(3) :: rdiff real(real_8), dimension(3) :: origin, xm real(real_8), dimension(3,3) :: box integer, dimension(3) :: n_points, n_points_r integer :: remainder integer :: i integer :: num_atoms character(len=*), parameter :: procedureN = 'mimic_ifc_init' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC allocate(wrapped_coordinates, source=tau) allocate(mimic_forces, mold=tau) mimic_forces = 0.0_real_8 origin(:) = 0.0_real_8 box(:,:) = 0.0_real_8 box(1,1) = cell_com%celldm(1) box(2,2) = cell_com%celldm(1) * cell_com%celldm(2) box(3,3) = cell_com%celldm(1) * cell_com%celldm(3) n_points(1) = spar%nr1s n_points(2) = spar%nr2s n_points(3) = spar%nr3s n_points_r(1) = fpar%kr1 n_points_r(2) = fpar%kr2s n_points_r(3) = fpar%kr3s call mimic_allocate_mm_struct(subsystems, & mimic_control%factors, & wrapped_coordinates, & mimic_forces, & ions0%na(1:mimic_control%num_quantum_species), & mimic_control%num_species, & mimic_control%covalent_radius, & sizes, & system_data) do i = 1, size(subsystems) CALL mp_bcast(subsystems(i)%num_atoms, parai%source, parai%cp_grp) end do if (paral%io_parent) then call mimic_init_constraints end if call mp_bcast(cotc0%mcnstr, & parai%io_source, & parai%cp_grp) if (cotc0%mcnstr > 0) then if (.not. paral%io_parent) then if (allocated(ntcnst)) then deallocate(ntcnst) deallocate(cnsval) deallocate(cnpar) deallocate(cnsval_dest) deallocate(grate) end if allocate(ntcnst(6, cotc0%mcnstr)) allocate(cnsval(cotc0%mcnstr)) allocate(cnpar(2, cotc0%mcnstr)) allocate(cnsval_dest(cotc0%mcnstr)) allocate(grate(cotc0%mcnstr)) end if call mp_bcast(ntcnst, & size(ntcnst), & parai%io_source, & parai%cp_grp) call mp_bcast(cnsval, & size(cnsval), & parai%io_source, & parai%cp_grp) end if if (paral%io_parent) then call mimic_collect_coordinates(subsystems) end if CALL mp_bcast(wrapped_coordinates, size(wrapped_coordinates), parai%io_source, parai%cp_grp) call mimic_allocate_qm_struct(quantum_fragment, ions0%na(1:mimic_control%num_quantum_species), & ions0%zv(1:mimic_control%num_quantum_species), & origin, box, n_points, n_points_r, rhoe, extf, & wrapped_coordinates, mimic_forces) call mimic_center_qm(rdiff, subsystems, quantum_fragment) xm(1) = box(1,1) * 0.5_real_8 xm(2) = box(2,2) * 0.5_real_8 xm(3) = box(3,3) * 0.5_real_8 call mimic_min_image(mimic_control%box, xm, subsystems) tau = wrapped_coordinates call initialize_tensors() if (mimic_control%long_range_coupling) then call mimic_compute_multipoles(quantum_fragment, mimic_control%multipole_order, & parap%NRXPL(parai%mepos,1), parap%NRXPL(parai%mepos,2)) call mp_sum(quantum_fragment%electronic_multipoles%values, & size(quantum_fragment%electronic_multipoles%values), & parai%allgrp) else call mimic_ifc_sort_fragments() end if #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_ifc_init !> re-center the QM subsystem and wrap the box with PBC subroutine mimic_update_coords(tau, c0, cm, nstates, npws, inyh) !> CPMD coordinate array real(real_8), intent(inout) :: tau(:,:,:) !> WF positions and velocities complex(kind=real_8), dimension(nkpt%ngwk, nstates, 2), intent(inout) :: c0, cm !> number of states and PWs integer, intent(in) :: nstates, npws !> Miller indices of PWs integer, dimension(:,:), intent(in) :: inyh real(real_8), dimension(3) :: grid_center real(real_8), dimension(3) :: rdiff real(real_8), dimension(3) :: rtrans, gdiff, xm real(real_8) :: gdot complex(real_8), dimension(:), allocatable :: eig_trans integer :: npw, nstate integer :: nk, nh real(real_8), dimension(3,3) :: box character(*), parameter :: procedureN = 'mimic_update_coords' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC grid_center(1) = spar%nr1s / 2 + 1 grid_center(2) = spar%nr2s / 2 + 1 grid_center(3) = spar%nr3s / 2 + 1 box(1,1) = cell_com%celldm(1) box(2,2) = cell_com%celldm(1) * cell_com%celldm(2) box(3,3) = cell_com%celldm(1) * cell_com%celldm(3) wrapped_coordinates = tau allocate(eig_trans(npws)) call mimic_center_qm(rdiff, subsystems, quantum_fragment) tau = wrapped_coordinates !\$OMP PARALLEL PRIVATE(npw, rtrans, gdiff, gdot, nstate) !\$OMP DO do npw = 1, npws rtrans = real(inyh(1:3,npw), real_8) - grid_center gdiff(1) = rtrans(1) * gvec_com%b1(1) & + rtrans(2) * gvec_com%b2(1) & + rtrans(3) * gvec_com%b3(1) gdiff(2) = rtrans(1) * gvec_com%b1(2) & + rtrans(2) * gvec_com%b2(2) & + rtrans(3) * gvec_com%b3(2) gdiff(3) = rtrans(1) * gvec_com%b1(3) & + rtrans(2) * gvec_com%b2(3) & + rtrans(3) * gvec_com%b3(3) gdot = parm%tpiba * dot_product(gdiff, rdiff) eig_trans(npw) = cmplx(cos(gdot), -sin(gdot), real_8) end do !\$OMP END DO !\$OMP DO do nstate = 1, nstates do npw = 1, npws c0(npw,nstate,1) = c0(npw,nstate,1) * eig_trans(npw) cm(npw,nstate,1) = cm(npw,nstate,1) * eig_trans(npw) end do !npw end do !nstate !\$OMP END DO if (cntl%tmdbo .and. cntl%textrap) then do nh = 1, cnti%mextra do nk = 1, nkpt%nkpnt !\$OMP DO do nstate = 1, nstates do npw = 1, npws cold(npw,nstate,nk,nh) = cold(npw,nstate,nk,nh) * eig_trans(npw) enddo enddo !\$OMP END DO enddo enddo endif !\$OMP END PARALLEL xm(1) = box(1,1) * 0.5_real_8 xm(2) = box(2,2) * 0.5_real_8 xm(3) = box(3,3) * 0.5_real_8 call mimic_min_image(mimic_control%box, xm, subsystems) clsaabox%mm_c_trans = clsaabox%mm_c_trans + rdiff #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_update_coords !> send atomic coordinates to clients subroutine mimic_ifc_send_coordinates real(real_8), dimension(3) :: shift character(*), parameter :: procedureN = 'mimic_send_coordinates' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC shift(1) = mimic_control%box(1,1) - cell_com%celldm(1) shift(2) = mimic_control%box(2,2) - cell_com%celldm(1) * cell_com%celldm(2) shift(3) = mimic_control%box(3,3) - cell_com%celldm(1) * cell_com%celldm(3) shift = -shift * 0.5_real_8 + clsaabox%mm_c_trans call mimic_translate(quantum_fragment, subsystems, -shift) call mimic_send_coords(subsystems) call mimic_translate(quantum_fragment, subsystems, shift) #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_ifc_send_coordinates !> determine QM and MM degrees of freedom subroutine mimic_subsystem_dof() integer :: i, j, ia, iat, is integer :: num_dof real(real_8) :: qm_constraints real(real_8) :: mm_constraints real(real_8) :: const character(*), parameter :: procedureN = 'mimic_subsystem_dof' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC ! compute constraints contributions to DoFs. qm_constraints = 0.0_real_8 mm_constraints = 0.0_real_8 if (cotc0%mcnstr > 0) then do i = 1, cotc0%mcnstr ! set prefactor depending on constraint type. if (ntcnst(1,i) == 1) then const = 0.5_real_8 else if (ntcnst(1,i) == 2) then const = 1.0_real_8/3.0_real_8 else if (ntcnst(1,i) == 3) then const = 0.25_real_8 else if (ntcnst(1,i) == 4) then const = 0.5_real_8 else if (ntcnst(1,i) == 5) then const = 0.25_real_8 else if (ntcnst(1,i) == 6) then const = 1.0_real_8 else if (ntcnst(1,i) == 7) then const = 1.0_real_8/3.0_real_8 else if (ntcnst(1,i) == 8) then const = 1.0_real_8 end if ! Loop over list of atoms in constraint ! Unused entries are supposed to be 0 do j = 2, 6 iat = ntcnst(j,i) if (iat <= 0 .or. iat > mimic_control%num_atoms) then const = 0.0_real_8 else if (iat <= mimic_control%num_quantum_atoms) then qm_constraints = qm_constraints + const else mm_constraints = mm_constraints + const endif end do end do end if num_dof = 0 do iat = 1, mimic_control%num_quantum_atoms ia = iatpt(1,i) is = iatpt(2,i) num_dof = num_dof + lskcor(1,iat) + lskcor(2,iat) + lskcor(3,iat) end do mimic_control%qm_dof = real(num_dof, kind=real_8) - qm_constraints num_dof = 0 do iat = mimic_control%num_quantum_atoms + 1, mimic_control%num_atoms ia = iatpt(1,iat) is = iatpt(2,iat) num_dof = num_dof + lskcor(1,iat) + lskcor(2,iat) + lskcor(3,iat) end do mimic_control%mm_dof = real(num_dof, kind=real_8) - mm_constraints if (paral%io_parent) then write(6,'(1x,a,t58,i8)') 'DEGREES OF FREEDOM FOR QM:', nint(mimic_control%qm_dof) write(6,'(1x,a,t58,i8)') 'DEGREES OF FREEDOM FOR MM:', nint(mimic_control%mm_dof) end if #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_subsystem_dof !> temperature of QM and MM subsystems separately (adapted from mm_localt) subroutine mimic_subsystem_temperatures(velocities) !> velocities of QM and MM subsystems real(real_8), dimension(:,:,:) :: velocities integer :: iat, ia, is real(real_8) :: qm_kinetic_energy real(real_8) :: mm_kinetic_energy real(real_8) :: qm_temperature real(real_8) :: mm_temperature character(*), parameter :: procedureN = 'mimic_subsystem_temperatures' integer :: isub call tiset(procedureN, isub) #ifdef __MIMIC qm_kinetic_energy = 0.0_real_8 do iat = 1, mimic_control%num_quantum_atoms ia = iatpt(1,iat) is = iatpt(2,iat) qm_kinetic_energy = qm_kinetic_energy + 0.5_real_8 * rmass%pma(is) * & (velocities(1,ia,is) * velocities(1,ia,is) + & velocities(2,ia,is) * velocities(2,ia,is) + & velocities(3,ia,is) * velocities(3,ia,is)) end do if (mimic_control%qm_dof > 0.1_real_8) then qm_temperature = factem * 2.0_real_8 * qm_kinetic_energy / mimic_control%qm_dof else qm_temperature = 0.0_real_8 end if mm_kinetic_energy = 0.0_real_8 do iat = mimic_control%num_quantum_atoms + 1, mimic_control%num_atoms ia = iatpt(1,iat) is = iatpt(2,iat) mm_kinetic_energy = mm_kinetic_energy + 0.5_real_8 * rmass%pma(is) * & (velocities(1,ia,is) * velocities(1,ia,is) + & velocities(2,ia,is) * velocities(2,ia,is) + & velocities(3,ia,is) * velocities(3,ia,is)) end do if (mimic_control%mm_dof > 0.1_real_8) then mm_temperature = factem * 2.0_real_8 * mm_kinetic_energy / mimic_control%mm_dof else mm_temperature = 0.0_real_8 endif if (paral%io_parent) then write(6,'(8x,a)') 'SUBSYSTEM TEMPERATURES:' write(6,'(8x,a,t58,f8.2)') 'T(QM) =', qm_temperature write(6,'(8x,a,t58,f8.2)') 'T(MM) =', mm_temperature end if #else call stopgm(procedureN, MIMIC_MISSING_ERROR, __LINE__, __FILE__) #endif call tihalt(procedureN, isub) end subroutine mimic_subsystem_temperatures end module mimic_wrapper . EOF cat << EOF > patch.24 diff ./src/mm_init_utils.mod.F90 ./src/mm_init_utils.mod.F90 458a IF (cntl%mimic) THEN CALL mimic_revert_dim() ENDIF . 393,394c IF (cntl%mimic) THEN mmdim%natq=mimic_control%num_quantum_atoms mmdim%naxq=mimic_control%max_quantum_atoms ELSE mmdim%natq=mmdim%natm mmdim%naxq=maxsys%nax END IF . 371a IF (cntl%mimic) THEN CALL mimic_save_dim() CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 34c USE system, ONLY: cntl,& maxsys . 14a USE mimic_wrapper, ONLY: mimic_save_dim,& mimic_switch_dim,& mimic_revert_dim,& mimic_control . EOF cat << EOF > patch.25 diff ./src/mp_interface.mod.F90 ./src/mp_interface.mod.F90 305a #ifdef __MIMIC call mimic_init(mp_comm_world) #endif . 15a #ifdef __MIMIC USE mimic_main, ONLY: mimic_init #endif . EOF cat << EOF > patch.26 diff ./src/mts_utils.mod.F90 ./src/mts_utils.mod.F90 96a ! Stops if MTS is used with MiMiC if (cntl%use_mts .and. cntl%mimic) then CALL stopgm(procedureN,'USE_MTS cannot be used with MIMIC',& __LINE__,__FILE__) end if . EOF cat << EOF > patch.27 diff ./src/nosalloc_utils.mod.F90 ./src/nosalloc_utils.mod.F90 287a IF (cntl%mimic) THEN CALL mimic_revert_dim() END IF . 60a IF (cntl%mimic) THEN CALL mimic_save_dim() CALL mimic_switch_dim(go_qm=.FALSE.) END IF . 10a USE mimic_wrapper, ONLY: mimic_save_dim,& mimic_revert_dim,& mimic_switch_dim . EOF cat << EOF > patch.28 diff ./src/nosepa_utils.mod.F90 ./src/nosepa_utils.mod.F90 898c CALL stopgm(procedureN,'NOSE + CONSTRAINTS',& . 871c ELSE DO is=1,mmdim%nspm DO ia=1,NAm(is) iat=iat+1 n1(iat+1)=lcttab(is,ia) ENDDO ENDDO END IF . 866,869c IF (cntl%mimic) THEN DO is=1,ions1%nsp DO ia=1,ions0%na(is) iat=iat+1 n1(iat+1)=lcttab(is,ia) ENDDO . 860c IF (cntl%mimic) THEN ALLOCATE(n1(3*maxsys%nax*ions1%nsp+1),STAT=ierr) ELSE ALLOCATE(n1(3*maxsys%nax*mmdim%nspm+1),STAT=ierr) END IF . 829,840c ! DO idxlct=1,loct%nloct ! ! WRITE(6,*)PARENT,'ATOMS IN THERMOSTATNO.',IDXLCT,NLCTMBM(IDXLCT) ! DO k=1,nlctmbm(idxlct) ! is=lctmemb(1,k,idxlct) ! ia=lctmemb(2,k,idxlct) ! iag=gratom((is-1)*maxsys%nax+ia) ! ! WRITE(6,'(I5)') IAG ! ENDDO ! ! PRINT * ! ENDDO . 771,773c IF (cntl%mimic) THEN ias=iatpt(1,ia) isp=iatpt(2,ia) ELSE isp=cpsp(ia) ias=cpat(ia) END IF . 170c IF (isos1%tisos.AND.(.NOT.lqmmm%qmmm.AND..NOT.cntl%mimic)) THEN . 37a iatpt,& . EOF cat << EOF > patch.29 diff ./src/pbc_utils.mod.F90 ./src/pbc_utils.mod.F90 109c IF (isos1%tclust.AND.mm_stat.AND..NOT.cntl%mimic) THEN . 107a IF (cntl%mimic) THEN a(:) = 0.0_real_8 a(1) = mimic_control%box(1,1) * 0.5_real_8 a(2) = mimic_control%box(2,2) * 0.5_real_8 a(3) = mimic_control%box(3,3) * 0.5_real_8 END IF . 7a USE mimic_wrapper, ONLY: mimic_control . EOF cat << EOF > patch.30 diff ./src/pcgrad_driver.mod.F90 ./src/pcgrad_driver.mod.F90 543a IF (cntl%mimic) THEN mimic_energy%qm_energy = ener_com%etot ener_com%etot = ener_com%etot & + ener_com%eext & + mimic_energy%qmmm_energy & + mimic_energy%mm_energy END IF . 228a IF (cntl%mimic) THEN CALL mimic_revert_dim() ENDIF . 117a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 113a IF (cntl%mimic) THEN CALL mimic_save_dim() CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 14a USE mimic_wrapper, ONLY: mimic_save_dim,& mimic_switch_dim,& mimic_revert_dim,& mimic_energy . EOF cat << EOF > patch.31 diff ./src/potfor_utils.mod.F90 ./src/potfor_utils.mod.F90 197a fion(1:3, 1:maxval(ions0%na(1:ions1%nsp)), 1:ions1%nsp) = fion_temp(:,:,:) . 190,192c fion_temp(1,ia,is)=fion_temp(1,ia,is)+REAL(ei123*txx)*omtp fion_temp(2,ia,is)=fion_temp(2,ia,is)+REAL(ei123*tyy)*omtp fion_temp(3,ia,is)=fion_temp(3,ia,is)+REAL(ei123*tzz)*omtp . 168c !\$omp reduction(+:FION_TEMP) . 149,151c fion_temp(1,ia,is)=fion_temp(1,ia,is)+REAL(ei123*txx)*omtp fion_temp(2,ia,is)=fion_temp(2,ia,is)+REAL(ei123*tyy)*omtp fion_temp(3,ia,is)=fion_temp(3,ia,is)+REAL(ei123*tzz)*omtp . 128c !\$omp reduction(+:FION_TEMP) . 107,109c fion_temp(1,ia,is)=fion_temp(1,ia,is)+REAL(ei123*gx)*omtp fion_temp(2,ia,is)=fion_temp(2,ia,is)+REAL(ei123*gy)*omtp fion_temp(3,ia,is)=fion_temp(3,ia,is)+REAL(ei123*gz)*omtp . 80,82c fion_temp(1,ia,is)=fion_temp(1,ia,is)+REAL(ei123*gx)*omtp fion_temp(2,ia,is)=fion_temp(2,ia,is)+REAL(ei123*gy)*omtp fion_temp(3,ia,is)=fion_temp(3,ia,is)+REAL(ei123*gz)*omtp . 57a allocate(fion_temp(3, maxval(ions0%na(1:ions1%nsp)), ions1%nsp)) fion_temp(:,:,:) = fion(1:3, 1:maxval(ions0%na(1:ions1%nsp)), 1:ions1%nsp) . 52a REAL(real_8), allocatable :: fion_temp(:,:,:) . EOF cat << EOF > patch.32 diff ./src/printp_utils.mod.F90 ./src/printp_utils.mod.F90 283c IF ((cotc0%mcnstr.GT.0.OR.cotr007%mrestr.GT.0).AND..NOT.cntl%new_constraints) THEN . 180c IF ((cotc0%mcnstr.GT.0.OR.cotr007%mrestr.GT.0).AND..NOT.cntl%new_constraints) THEN . EOF cat << EOF > patch.33 diff ./src/proppt_utils.mod.F90 ./src/proppt_utils.mod.F90 1026a IF (cntl%mimic) THEN CALL mimic_revert_dim() ENDIF . 1013a IF (paral%io_parent.AND.cntl%mimic) THEN CALL mimic_ifc_send_coordinates() CALL mimic_ifc_collect_energies() CALL mimic_ifc_collect_forces() END IF . 867a IF (cntl%mimic) THEN CALL mimic_ifc_init(rhoe, extf, tau0) END IF . 803a IF (cntl%mimic) THEN CALL mimic_ifc_init(rhoe, extf, tau0) END IF . 681a IF (cntl%mimic) THEN CALL mimic_ifc_init(rhoe, extf, tau0) END IF . 574a IF (cntl%mimic) THEN CALL mimic_ifc_init(rhoe, extf, tau0) END IF . 546a IF (cntl%mimic) THEN CALL mimic_ifc_init(rhoe, extf, tau0) END IF . 529a IF (cntl%mimic) THEN CALL mimic_ifc_init(rhoe, extf, tau0) END IF . 510a IF (cntl%mimic) THEN CALL mimic_ifc_init(rhoe, extf, tau0) END IF . 497a IF (cntl%mimic) THEN CALL mimic_ifc_init(rhoe, extf, tau0) END IF . 484a IF (cntl%mimic) THEN CALL mimic_ifc_init(rhoe, extf, tau0) END IF . 461a IF (cntl%mimic) THEN ALLOCATE(extf(fpar%kr1*fpar%kr2s*fpar%kr3s),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) CALL mimic_ifc_init(rhoe, extf, tau0) mimic_control%update_potential = .TRUE. END IF . 436a IF (cntl%mimic) THEN CALL mimic_ifc_init(rhoe, extf, tau0) END IF . 346a IF (cntl%mimic) THEN CALL mimic_ifc_init(rhoe, extf, tau0) END IF . 251a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 247a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 225a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 222,224c IF (.NOT.cntl%mimic) THEN ALLOCATE(taup(3,maxsys%nax,maxsys%nsx),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) ENDIF . 179a IF (cntl%mimic) THEN CALL mimic_save_dim() CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 69a USE mimic_wrapper, ONLY: mimic_ifc_collect_energies,& mimic_ifc_collect_forces,& mimic_ifc_sort_fragments,& mimic_ifc_init,& mimic_revert_dim,& mimic_save_dim,& mimic_ifc_send_coordinates,& mimic_sum_forces,& mimic_switch_dim,& mimic_control . 25a USE efld, ONLY: extf . EOF cat << EOF > patch.34 diff ./src/ratom_utils.mod.F90 ./src/ratom_utils.mod.F90 749a ENDIF IF (cntl%mimic) THEN maxsys%nsx = mimic_control%num_quantum_species maxsys%nax = mimic_control%max_quantum_atoms . 395a IF (cntl%mimic) maxsys%nsx = mimic_control%num_species . 387a IF (cntl%mimic) maxsys%nsx = mimic_control%num_quantum_species . 161a IF (cntl%mimic) THEN NSX_q = mimic_control%num_quantum_species maxsys%nsx = mimic_control%num_species maxsys%nax = mimic_control%max_atoms ENDIF . 66c USE system, ONLY: cntl,& maxsp,& . 38a USE mimic_wrapper, ONLY: mimic_control . EOF cat << EOF > patch.35 diff ./src/rv30_utils.mod.F90 ./src/rv30_utils.mod.F90 194a ENDIF IF (cntl%mimic) THEN CALL mimic_revert_dim() . 123a IF (cntl%mimic) THEN CALL mimic_save_dim() CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 41a USE mimic_wrapper, ONLY: mimic_save_dim,& mimic_switch_dim,& mimic_revert_dim . EOF cat << EOF > patch.36 diff ./src/rwfopt_utils.mod.F90 ./src/rwfopt_utils.mod.F90 973a IF (cntl%mimic) THEN CALL mimic_revert_dim() ENDIF . 861a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 853a IF (cntl%mimic .AND. .NOT. mimic_control%tot_scf_energy) THEN IF (paral%io_parent) THEN CALL mimic_ifc_collect_energies() END IF END IF . 852a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 847a IF (cntl%mimic) THEN CALL mimic_sum_forces(fion) ENDIF . 830a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 665a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 662a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 523a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 501a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 320a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 319a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 312a IF (cntl%mimic) THEN CALL mimic_update_coords(tau0, c0, dummy, crge%n, ncpw%ngw, inyh) IF (paral%io_parent) THEN CALL mimic_ifc_send_coordinates() IF (mimic_control%tot_scf_energy) THEN CALL mimic_ifc_collect_energies() END IF END IF END IF . 250a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 195c IF (.NOT.cntl%bsymm.AND..NOT.cntl%mimic) THEN . 193a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 192a IF (cntl%mimic) THEN ALLOCATE(dummy, mold=c0) ALLOCATE(extf(fpar%kr1*fpar%kr2s*fpar%kr3s),STAT=ierr) IF (ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) CALL zeroing(extf) !,kr1*kr2s*kr3s) CALL mimic_ifc_init(rhoe, extf, tau0) CALL zeroing(taup) !,3*maxsys%nax*maxsys%nsx) CALL zeroing(fion) !,3*maxsys%nax*maxsys%nsx) mimic_control%update_potential = .TRUE. END IF . 153a IF (cntl%mimic) THEN CALL mimic_save_dim() ENDIF . 150a COMPLEX(real_8), DIMENSION(:,:,:), ALLOCATABLE :: dummy . 81c cnti, cntl, cntr, fpar, locpot2, maxsys, nacc, ncpw, nkpt, parm, spar, parap . 58c USE mp_interface, ONLY: mp_sum,& mp_bcast . 53a USE mimic_wrapper, ONLY: mimic_ifc_collect_energies,& mimic_ifc_collect_forces,& mimic_ifc_sort_fragments,& mimic_ifc_init,& mimic_revert_dim,& mimic_save_dim,& mimic_ifc_send_coordinates,& mimic_sum_forces,& mimic_switch_dim,& mimic_update_coords,& mimic_control . 22c USE efld, ONLY: extf,& textfld . 18a USE cppt, ONLY: inyh . EOF cat << EOF > patch.37 diff ./src/setirec_utils.mod.F90 ./src/setirec_utils.mod.F90 193c irec(irec_ctrans) = isetone(cntl%tqmmm.OR.cntl%mimic) . 124c irec(irec_ctrans) = isetone(cntl%tqmmm.OR.cntl%mimic) . EOF cat << EOF > patch.38 diff ./src/setsys_utils.mod.F90 ./src/setsys_utils.mod.F90 1622a IF (cntl%mimic) THEN CALL mimic_revert_dim() ENDIF . 427a IF (cntl%mimic) THEN CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 358c .AND..NOT.cntl%tmdeh.AND..NOT.cntl%tpspec& .AND..NOT.cntl%mimic) THEN . 161a IF (cntl%mimic) THEN CALL mimic_save_dim() CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 56a USE mimic_wrapper, ONLY: mimic_save_dim,& mimic_switch_dim,& mimic_revert_dim . EOF cat << EOF > patch.39 diff ./src/shake_utils.mod.F90 ./src/shake_utils.mod.F90 30a SUBROUTINE init_constraints(ntcnst, cval, na, nsp) ! CPMD constraint map integer, dimension(:,:), intent(in) :: ntcnst ! Reference values for constraints real(real_8), dimension(:), intent(in) :: cval ! Number of atoms per species integer, dimension(:), intent(in) :: na ! Number of atomic species integer, intent(in) :: nsp character(*), PARAMETER :: procedureN = 'init_constraints' integer :: isub call tiset(procedureN,isub) const_matrix = build_constraints(ntcnst, cnsval, ions0%na, ions1%nsp) xlagr(:) = 0.0_real_8 ylagr(:) = 0.0_real_8 call tihalt(procedureN,isub) END SUBROUTINE init_constraints SUBROUTINE do_shake(tau0, taup, velp) !> CPMD coordinate matrix real(real_8), dimension(:,:,:) :: tau0 real(real_8), dimension(:,:,:) :: taup !> CPMD velocity matrix real(real_8), dimension(:,:,:) :: velp character(*), PARAMETER :: procedureN = 'do_shake' integer :: isub call tiset(procedureN,isub) xlagr(:) = ylagr(:) call new_shake(const_matrix, xlagr, tau0, taup, velp, dt_ions, dtb2mi) call tihalt(procedureN,isub) END SUBROUTINE do_shake SUBROUTINE do_rattle(tau0, velp) !> CPMD coordinate matrix real(real_8), dimension(:,:,:) :: tau0 !> CPMD velocity matrix real(real_8), dimension(:,:,:) :: velp character(*), PARAMETER :: procedureN = 'do_rattle' integer :: isub call tiset(procedureN,isub) ylagr(:) = 0.0_real_8 call new_rattle(const_matrix, ylagr, tau0, velp, dt_ions, dtb2mi) call tihalt(procedureN,isub) END SUBROUTINE do_rattle . 26a PUBLIC :: init_constraints PUBLIC :: do_shake PUBLIC :: do_rattle . 25a type(constraint_matrix) :: const_matrix . 5a USE constraint USE constraint_utils . EOF cat << EOF > patch.40 diff ./src/SOURCES ./src/SOURCES 272c dftd3_driver.mod.F90 mimic_wrapper.mod.F90 . 140,141c setbsstate_utils.mod.F90 io_utils.mod.F90 atoms_utils.mod.F90 dynit_utils.mod.F90 constraint.mod.F90 constraint_utils.mod.F90 \\ shake_utils.mod.F90 rattle_utils.mod.F90 resetac_utils.mod.F90 dispp_utils.mod.F90 \\ . EOF cat << EOF > patch.41 diff ./src/sysin_utils.mod.F90 ./src/sysin_utils.mod.F90 1447a IF (cntl%mimic.AND.parm%ibrav/=0) THEN WRITE(output_unit,'(A)') ' ERROR: SYMMETRY IBRAV =',parm%ibrav,' NOT ALLOWED IN MiMiC' WRITE(output_unit,'(A)') ' A VALUE OF 0 (ISOLATED) IS REQUIRED' CALL stopgm(procedureN, 'Cell dimensions incompatible with MiMiC', & __LINE__, __FILE__) ENDIF . EOF cat << EOF > patch.42 diff ./src/system.mod.F90 ./src/system.mod.F90 792a REAL(real_8), DIMENSION(2) :: anneal_factors = HUGE(0.0_real_8) . 675a INTEGER :: shake_maxstep = HUGE(0) INTEGER :: shake_cg_iter = HUGE(0) . 565a LOGICAL :: mimic = .FALSE. LOGICAL :: new_constraints = .FALSE. LOGICAL :: pbicgstab = .FALSE. LOGICAL :: anneal_dual = .FALSE. . EOF cat << EOF > patch.43 diff ./src/updrho_utils.mod.F90 ./src/updrho_utils.mod.F90 658c IF (cntl%mimic) THEN mimic_energy%qm_energy = ener_com%etot ener_com%etot = ener_com%etot & + ener_com%eext & + mimic_energy%qmmm_energy & + mimic_energy%mm_energy END IF ! Check total energy difference . 55a USE mimic_wrapper, ONLY: mimic_energy . EOF cat << EOF > patch.44 diff ./src/wrener_utils.mod.F90 ./src/wrener_utils.mod.F90 685a ELSE IF (cntl%mimic) THEN IF (paral%io_parent) THEN defenergy_mimic = defenergy(1:len(trim(adjustl(defenergy)))-1)//'+MM+QM/MM)' WRITE(6,'(/,1X,A,T25,A,T41,F20.8,A)')& defenergy_mimic,'TOTAL ENERGY =',ener_com%etot,' A.U.' WRITE(6,'(1X,A,T22,A,T41,F20.8,A)')& defenergy, 'TOTAL QM ENERGY =', mimic_energy%qm_energy, ' A.U.' WRITE(6,'(1X,A,T22,A,T41,F20.8,A)')& '(MM)', 'TOTAL MM ENERGY =', mimic_energy%mm_energy, ' A.U.' WRITE(6,'(1X,A,T19,A,T41,F20.8,A)')& '(QM/MM)', 'TOTAL QM/MM ENERGY =', mimic_energy%qmmm_energy + ener_com%eext, ' A.U.' ENDIF . 541a CHARACTER(len=25) :: defenergy_mimic . 332c IF (.NOT.cntl%new_constraints) THEN CALL cnstpr END IF . 27a USE mimic_wrapper, ONLY: mimic_control,& mimic_energy . EOF cat << EOF > patch.45 diff ./src/wrgeo_utils.mod.F90 ./src/wrgeo_utils.mod.F90 205a ENDIF IF (cntl%mimic) THEN CALL mimic_revert_dim() . 163a IF (cntl%mimic) THEN CALL mimic_save_dim() CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 122a IF (cntl%mimic) THEN CALL mimic_revert_dim() ENDIF . 91a IF (cntl%mimic) THEN CALL mimic_revert_dim() ENDIF . 89a IF (cntl%mimic) THEN CALL mimic_save_dim() CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 76a IF (cntl%mimic) THEN CALL mimic_revert_dim() ENDIF . 45a IF (cntl%mimic) THEN CALL mimic_save_dim() CALL mimic_switch_dim(go_qm=.TRUE.) ENDIF . 24c USE system, ONLY: cnti,& cntl . 11a USE mimic_wrapper, ONLY: mimic_save_dim,& mimic_switch_dim,& mimic_revert_dim . EOF cat << EOF > patch.46 diff ./src/wv30_utils.mod.F90 ./src/wv30_utils.mod.F90 227a ENDIF IF (cntl%mimic) THEN CALL mimic_revert_dim() . 127a IF (cntl%mimic) THEN CALL mimic_save_dim() CALL mimic_switch_dim(go_qm=.FALSE.) ENDIF . 44a USE mimic_wrapper, ONLY: mimic_save_dim,& mimic_switch_dim,& mimic_revert_dim . EOF # Progress bar from https://www.baeldung.com/linux/command-line-progress-bar bar_size=40 bar_char_done="#" bar_char_todo="-" bar_percentage_scale=2 function show_progress { current="$1" total="$2" # calculate the progress in percentage percent=$(bc <<< "scale=$bar_percentage_scale; 100 * $current / $total" ) # The number of done and todo characters done=$(bc <<< "scale=0; $bar_size * $percent / 100" ) todo=$(bc <<< "scale=0; $bar_size - $done" ) # build the done and todo sub-bars done_sub_bar=$(printf "%${done}s" | tr " " "${bar_char_done}") todo_sub_bar=$(printf "%${todo}s" | tr " " "${bar_char_todo}") # output the bar echo -ne "\rProgress : [${done_sub_bar}${todo_sub_bar}] ${percent}%" if [ $total -eq $current ]; then echo -e "\nDONE" fi } echo "Patching CPMD" echo "Adding MiMiC interface to CPMD" for i in `seq -w 0 46` do filename=`head -n 1 patch.$i | awk '{print $2}'` patch -e -p0 $filename < patch.$i rm patch.$i show_progress $i 46 done echo "MiMiC patch has been applied"