diff --git a/configure/MIMIC-LINUX-GCC b/configure/MIMIC-LINUX-GCC new file mode 100644 index 0000000..d7a5851 --- /dev/null +++ b/configure/MIMIC-LINUX-GCC @@ -0,0 +1,23 @@ + 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 " diff --git a/configure/MIMIC-LINUX-INTEL b/configure/MIMIC-LINUX-INTEL new file mode 100644 index 0000000..72e77ea --- /dev/null +++ b/configure/MIMIC-LINUX-INTEL @@ -0,0 +1,25 @@ + 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 " diff --git a/doc/manual/manual.tex b/doc/manual/manual.tex index d049e1f..44df685 100644 --- a/doc/manual/manual.tex +++ b/doc/manual/manual.tex @@ -657,6 +657,8 @@ \subsubsection{Version 4.3} \subsubsection{Version 4.5} The Minnesota non-separable exchange functional families (MN12, MN15) as well as the revM06-L functional are now available. +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. % XXX tutorial added. \clearpage @@ -1117,6 +1119,10 @@ \subsection{Input Sections}\label{sections} \\ \> \&MTS ... \> \&END \> $\leftrightarrow$ \> Parameters for the Multiple Time-Step MD scheme. \\ +\\ +\> \&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.\\ \end{tabbing} A detailed discussion of the different keywords will be given in the @@ -1133,7 +1139,7 @@ \subsection{List of Keywords by Sections} \refkeyword{ALEXANDER MIXING}\options{}{}{} \refkeyword{ALLTOALL}\options{\{SINGLE,DOUBLE\}}{}{} \refkeyword{ANDERSON MIXING}\options{}{}{} -\refkeyword{ANNEALING}\options{\{IONS,ELECTRONS,CELL\}}{}{} +\refkeyword{ANNEALING}\options{\{IONS,ELECTRONS,CELL\}}{[DUAL]}{} \refkeyword{BENCHMARK}\options{}{}{} \refkeyword{BERENDSEN}\options{\{IONS,ELECTRONS,CELL\}}{}{} \refkeyword{BFGS}\options{}{}{} @@ -1165,6 +1171,7 @@ \subsection{List of Keywords by Sections} \refkeyword{DEBUG NOACC}\options{}{}{} \refkeyword{DIIS MIXING}\options{}{}{} \refkeyword{DIPOLE DYNAMICS}\options{\{SAMPLE,WANNIER\}}{}{} +\refkeyword{DISABLE CONSTRAINTS}\options{}{}{} \refkeyword{DISTRIBUTED LINALG}\options{\{ON,OFF\}}{}{} \refkeyword{DISTRIBUTE FNL}\options{}{}{} \refkeyword{ELECTRONIC SPECTRA}\options{}{}{} @@ -1207,6 +1214,7 @@ \subsection{List of Keywords by Sections} \refkeyword{MAXITER}\options{}{}{} \refkeyword{MAXSTEP}\options{}{}{} \refkeyword{MEMORY}\options{}{[SMALL, BIG]}{} +\refkeyword{MIMIC}\options{}{}{} \refkeyword{MIRROR}\options{}{}{} \refkeyword{MIXDIIS}\options{}{}{} \refkeyword{MIXSD}\options{}{}{} @@ -1214,6 +1222,7 @@ \subsection{List of Keywords by Sections} \refkeyword{MOLECULAR DYNAMICS}\options{\{CP, BO, EH, PT, CLASSICAL, FILE [XYZ, NSKIP=N, NSAMPLE=M]\}}{}{} \refkeyword{MOVERHO}\options{}{}{} \refkeyword{MOVIE}\options{}{[OFF, SAMPLE]}{} +\refkeyword{NEW CONSTRAINTS}\options{}{[PBICGSTAB]}{} \refkeyword{NOGEOCHECK}\options{}{}{} \refkeyword{NONORTHOGONAL ORBITALS}\options{}{[OFF]}{} \refkeyword{NOSE}\options{\{IONS, ELECTRONS, CELL\}}{[ULTRA,MASSIVE,CAFES]}{} @@ -1253,6 +1262,8 @@ \subsection{List of Keywords by Sections} \refkeyword{RHOOUT}\options{}{[BANDS,SAMPLE=nrhoout]}{} \refkeyword{ROKS}\options{}{\{SINGLET, TRIPLET\},\{DELOCALIZED, LOCALIZED, GOEDECKER\}}{} \refkeyword{SCALED MASSES}\options{}{[OFF]}{} +\refkeyword{SHAKE\_CG\_ITER}\options{}{}{} +\refkeyword{SHAKE\_MAXSTEP}\options{}{}{} \refkeyword{SHIFT POTENTIAL}\options{}{}{} \refkeyword{SOC}\options{}{}{} \refkeyword{SPLINE}\options{}{[POINTS, QFUNCTION, INIT, RANGE]}{} @@ -1660,6 +1671,22 @@ \subsection{List of Keywords by Sections} \refspekeyword{LOW\_LEVEL\_FORCES}{LOW LEVEL FORCES}\options{\{DFT, EXTERNAL\}}{}{} \refspekeyword{HIGH\_LEVEL\_FORCES}{HIGH LEVEL FORCES}\options{\{DFT, EXTERNAL\}}{}{} % +\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{}{}{} +% \clearpage % %--------------------------------------------------------------------- @@ -1757,9 +1784,14 @@ \subsection{Alphabetic List of Keywords} Not supported for \refkeyword{QMMM} calculations. } -\keyword{ANNEALING}{\{IONS,ELECTRONS,CELL\}}{}{}{\&CPMD} +\keyword{ANNEALING}{\{IONS,ELECTRONS,CELL\}}{DUAL}{}{\&CPMD} \desc{Scale the ionic, electronic, or cell velocities every - time step. The scaling factor is read from the next line.} + 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.} \keyword{ATOMIC CHARGES}{}{}{}{\&ATOMS} \desc{Changes the default charge (0) of the atoms for the initial guess to @@ -1928,6 +1960,11 @@ \subsection{Alphabetic List of Keywords} % The keyword has to appear after \refkeyword{FREE ENERGY FUNCTIONAL}.} +\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.} + \keyword{BOX WALLS}{}{}{}{\&CPMD} \desc{The thickness parameter for soft, reflecting QM-box walls is read from the next line. This keyword allows to reverse the @@ -2423,6 +2460,13 @@ \subsection{Alphabetic List of Keywords} $|g + k|^2 < E_{cut}$ instead of $|g|^2 < E_{cut}$. This is the default.} +\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.} + \keyword{CZONES}{}{[SET]}{}{\&CPMD} \desc{Activates convergence zones for the wavefunction during the \refkeyword{CDFT} constraint minimisation. If SET is set the @@ -2606,6 +2650,11 @@ \subsection{Alphabetic List of Keywords} % {\bf Default} is to use the real-space algorithm.} +\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.} + \keyword{DISCARD}{}{[OFF, PARTIAL, TOTAL, LINEAR]}{}{\&RESP} \desc{Request to discard trivial modes in vibrational analysis from linear response (both \refkeyword{PHONON} and \refkeyword{LANCZOS}).\\ @@ -2978,6 +3027,18 @@ \subsection{Alphabetic List of Keywords} \desc{The state for which the forces are calculated is read from the next line. Default is for state 1.} +\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.} + \keyword{FREE ENERGY FUNCTIONAL}{}{}{}{\&CPMD} \desc{Calculates the electronic free energy using free energy density functional~\cite{Alavi94,PSil,mbaops} @@ -3740,6 +3801,16 @@ \subsection{Alphabetic List of Keywords} These wavefunction can be printed with the keyword {\bf RHOOUT} specified in the section \&CPMD section.} +\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.} + \keyword{LR KERNEL}{}{}{functional}{\&DFT} \desc{Use another functional for the linear response kernel. To be used like \refkeyword{FUNCTIONAL}.} @@ -3839,6 +3910,11 @@ \subsection{Alphabetic List of Keywords} \desc{Initiate Metadynamics (see section~\ref{sec:meta} for more information on the available options and the input format).} +\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.} + \keyword{MIRROR}{}{}{}{\&CPMD} \desc{Write the input file to the output.} @@ -3940,6 +4016,11 @@ \subsection{Alphabetic List of Keywords} {\bf Default} is {\bf not} to write a movie file.} +\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.} + \keyword{MULTIPLICITY}{}{}{}{\&SYSTEM} \desc{This keyword only applies to LSD calculations.\\ The multiplicity (2$S$+1) is read from the next line.\\ @@ -3964,6 +4045,19 @@ \subsection{Alphabetic List of Keywords} \keyword{NEQUI}{}{}{}{\&PATH} \desc{Number of equilibration steps discarded to calculate the mean force.} +\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.} + \keyword{NEWCODE}{}{}{}{\&DFT}% \desc{This keyword will be deprecated in a future release. Both \textbf{NEWCODE} and \textbf{OLDCODE} have been replaced by the new @@ -4114,7 +4208,8 @@ \subsection{Alphabetic List of Keywords} \end{tabular} The parser for the atom ranges uses either the CPMD ordering or - the GROMOS ordering in case of classical or QM/MM runs. Multiple + 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 ranges may be specified for the same thermostat. Atoms belonging to the same CPMD constraint or the same solvent molecule in QM/MM runs must belong to the same local thermostat. @@ -4243,6 +4338,36 @@ \subsection{Alphabetic List of Keywords} \spekeyword{OUTPUT}{}{[ALL, GROUPS, PARENT]}{}{\&PATH}{OUTPUT PATH} \desc{Idem as above, here for Mean Free Energy Path runs.} +\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.} \spekeyword{PARA\_BUFF\_SIZE}{}{}{}{\&CPMD}{PARA BUFF SIZE} \desc{ Set the buffer size for parallel operation (sum, ...). @@ -4274,6 +4399,10 @@ \subsection{Alphabetic List of Keywords} With the additional keywork { \bf SHOCK} } a MD simulation using the multiscale shock method \cite{shock} is performed. +\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.} + \keyword{PATH INTEGRAL}{}{}{}{\&CPMD} \desc{Perform a {\bf path integral molecular dynamics} calculation~\cite{Marx94,Marx96}.\\ @@ -5052,6 +5181,14 @@ \subsection{Alphabetic List of Keywords} } +\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}).} + \keyword{SHIFT POTENTIAL}{}{}{}{\&CPMD} \desc{After this keyword, useful in Hamiltonian diagonalization, the shift value $V_{\rm shift}$ must be provided in the next line. @@ -5220,6 +5357,14 @@ \subsection{Alphabetic List of Keywords} to a set of atoms. The number of atoms and a list of the selected atoms has to be given on the next lines.} +\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.} + + \keyword{SUBTRACT}{}{[COMVEL, ROTVEL]}{}{\&CPMD} \desc{If COMVEL is selected, the total momentum of the system is removed, if ROTVEL is selected the global angular momentum of the @@ -5394,7 +5539,7 @@ \subsection{Alphabetic List of Keywords} \keyword{TEMPERATURE}{}{[RAMP]}{}{\&CPMD} \desc{The {\bf initial temperature} in Kelvin of the {\bf system} is read from the next line. With the additional keyword {\bf RAMP} the temperature can be linearly ramped to a target - value and two more numbers are read, the ramping target temperature in Kelvin and the + value and two more numbers are read, the ramping target temperature in Kelvin and the ramping speed in Kelvin per atomic time unit (to get the change per timestep you have to multiply it with the value of \refkeyword{TIMESTEP}). Note that this ramping affects the target temperatures for \refkeyword{TEMPCONTROL}, \refkeyword{BERENDSEN} and the global \refkeyword{NOSE} thermostats.} @@ -5453,6 +5598,14 @@ \subsection{Alphabetic List of Keywords} {\bf Default} is a factor of 1, i.e. the MD will be identical to a Velocity-Verlet MD with the high level forces. } +\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.} + \keyword{TRACE}{}{[ALL,MASTER]}{}{\&CPMD} \desc{ Activate the tracing of the procedures. {\sl ALL} specifies that all the mpi tasks are traced. {\sl ALL} specifies that only the master is traced. @@ -11985,6 +12138,18 @@ \section*{References} % R. Le Sar, J.~Phys.~Chem. {\bf 88}, 4272 (1984); % R. Le Sar, J.~Chem.~Phys. {\bf 86}, 1485 (1987). +\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} + \end{thebibliography} % %--------------------------------------------------------------------- diff --git a/scripts/configure.sh b/scripts/configure.sh index f4ba34a..4db95f8 100755 --- a/scripts/configure.sh +++ b/scripts/configure.sh @@ -57,7 +57,8 @@ Description of options: -verbose (-v) Display the name of files for the dependencies -coverage (-c) With GCC compiler only, allows for specific configuration files, to generate the code coverage/profiling during an execution - -qmmm Generates a makefile for QMMM + -mimic Enable interface to MiMiC (requires MPI and that MiMiC is available) + -qmmm Generates a makefile for QMMM -iphigenie Support for external interface to iphigenie -omp Enables the use of OMP instructions (if the config file allows that) OMP3 is in general triggered by the config keyword OMP3_DISABLED, @@ -130,6 +131,10 @@ do -disable_omp3) omp3=0 ;; + -mimic|-m) + mimic=1 + echo "** Enabling MiMiC interface (if MiMiC is available)" >&2 + ;; -qmmm|-q) qmmm=1 echo "** Enabling QM/MM (if Gromos and MM_Interface modules will be available)" >&2 @@ -241,6 +246,25 @@ fi CPPFLAGS_GROMOS='-DEWALD -DEWATCUT -DHAT_SHAPE -DUNPACKED_GRID' +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 #QM/MM compilation setup if [ $qmmm ]; then if [ -f ${MOD_DIR}/QMMM_SOURCES ]; then @@ -852,6 +876,11 @@ do SkipInclude["mpi"] = 0; SkipInclude["rhjsx.inc"] = 0; SkipInclude["uhjsx.inc"] = 0; + SkipInclude["mimic_precision"] = 0; + SkipInclude["mimic_constants"] = 0; + SkipInclude["mimic_types"] = 0; + SkipInclude["mimic_tensors"] = 0; + SkipInclude["mimic_main"] = 0; if (qmmm != 1) { SkipInclude["coordsz"] = 0; } diff --git a/src/SOURCES b/src/SOURCES index 3dfa594..863c4a5 100644 --- a/src/SOURCES +++ b/src/SOURCES @@ -137,8 +137,8 @@ SRC_MODULES = kinds.mod.F90 machine.mod.F90 timer.mod.F90 error_handling.mod.F90 meta_colvar_inp_utils.mod.F90 meta_exlagr_methods.mod.F90 meta_cv_utils.mod.F90 \ meta_cv_qmmm_utils.mod.F90 tst2min_utils.mod.F90 tst2min_inp_utils.mod.F90 chain_dr_utils.mod.F90 meta_cell_utils.mod.F90 \ meta_exl_mult_utils.mod.F90 meta_ex_mul_util_utils.mod.F90 meta_localizespin_utils.mod.F90 meta_colvar_util_utils.mod.F90 \ - setbsstate_utils.mod.F90 io_utils.mod.F90 atoms_utils.mod.F90 dynit_utils.mod.F90 shake_utils.mod.F90 \ - rattle_utils.mod.F90 resetac_utils.mod.F90 dispp_utils.mod.F90 \ + 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 \ noseinit_utils.mod.F90 nospinit_utils.mod.F90 noseng_utils.mod.F90 nosepa_utils.mod.F90 noseup_utils.mod.F90 \ enosmove_utils.mod.F90 pnosmove_utils.mod.F90 \ vdw_utils.mod.F90 chksym_utils.mod.F90 symtrz_utils.mod.F90 multtb_utils.mod.F90 molsym_utils.mod.F90 \ @@ -269,4 +269,4 @@ SRC_MODULES = kinds.mod.F90 machine.mod.F90 timer.mod.F90 error_handling.mod.F90 bicanonicalCpmd.mod.F90 bicanonicalConfig.mod.F90 bicanonicalInputReader.mod.F90 bicanonicalCalculationConfig.mod.F90 bicanonicalCalculation.mod.F90 \ pi_prpt_utils.mod.F90 pi_npt_cpmd_utils.mod.F90 pi_stress_utils.mod.F90 pi_npt_bomd_utils.mod.F90 \ nvml_interfaces.mod.F90 nvml_types.mod.F90 nvml_utils.mod.F90 mts_utils.mod.F90 interface_utils.mod.F90 \ - dftd3_driver.mod.F90 + dftd3_driver.mod.F90 mimic_wrapper.mod.F90 diff --git a/src/anneal_utils.mod.F90 b/src/anneal_utils.mod.F90 index 62e3acf..7e5a96e 100644 --- a/src/anneal_utils.mod.F90 +++ b/src/anneal_utils.mod.F90 @@ -5,6 +5,7 @@ MODULE anneal_utils s_ekinpp USE ions, ONLY: ions1 USE kinds, ONLY: real_8 + USE mimic_wrapper, ONLY: mimic_control USE nose, ONLY: glib USE parac, ONLY: paral USE pimd, ONLY: ipcurr,& @@ -51,19 +52,48 @@ SUBROUTINE anneal(velp,cm,nstate,htvel) ! == ANNEALING (IONS) Fixed rescaling factor (anner) == ! ==--------------------------------------------------------------== IF (cntl%annei.AND.(.NOT.skipanic)) THEN - alfap=cntr%anneri**(0.25_real_8) + IF (cntl%anneal_dual) THEN + alfap = cntr%anneal_factors(1)**(0.25_real_8) #if defined(__VECTOR) - !$omp parallel do private(I,IS,IA) + !$omp parallel do private(I,IS,IA) #else - !$omp parallel do private(I,IS,IA) schedule(static) + !$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 + 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 ENDIF ! ==--------------------------------------------------------------== ! == ANNEALING (ELECTRONS) Fixed rescaling factor (anner) == diff --git a/src/constraint.mod.F90 b/src/constraint.mod.F90 new file mode 100644 index 0000000..8db8afd --- /dev/null +++ b/src/constraint.mod.F90 @@ -0,0 +1,126 @@ +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 diff --git a/src/constraint_utils.mod.F90 b/src/constraint_utils.mod.F90 new file mode 100644 index 0000000..b0eee9b --- /dev/null +++ b/src/constraint_utils.mod.F90 @@ -0,0 +1,663 @@ +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 diff --git a/src/control_def_utils.mod.F90 b/src/control_def_utils.mod.F90 index fad44cb..13ef575 100644 --- a/src/control_def_utils.mod.F90 +++ b/src/control_def_utils.mod.F90 @@ -665,6 +665,13 @@ SUBROUTINE control_def iface1%intwrite = .FALSE. intfn = "interface.bin" ! ==--------------------------------------------------------------== + 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 RETURN END SUBROUTINE control_def ! ================================================================== diff --git a/src/control_pri_utils.mod.F90 b/src/control_pri_utils.mod.F90 index b48000f..95959e6 100644 --- a/src/control_pri_utils.mod.F90 +++ b/src/control_pri_utils.mod.F90 @@ -513,8 +513,15 @@ SUBROUTINE control_pri(unknown,iunk,fformat) 'EVERY ',cnti%imovie,' STEPS IN MOVIE FORMAT' ENDIF IF (cntl%annei) THEN - WRITE(6,'(A,T46,A8,F12.6)')& - ' SIMULATED ANNEALING OF IONS WITH ','ANNERI =',cntr%anneri + 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 ENDIF IF (cntl%annee) THEN WRITE(6,'(A,T46,A8,F12.6)')& diff --git a/src/control_utils.mod.F90 b/src/control_utils.mod.F90 index 4fa1b5e..27d99ad 100644 --- a/src/control_utils.mod.F90 +++ b/src/control_utils.mod.F90 @@ -3159,8 +3159,23 @@ SUBROUTINE control previous_line = line READ(iunit,'(A)',iostat=ierr) line IF ( keyword_contains(previous_line,'IONS') ) THEN - cntl%annei=.TRUE. - CALL readsr(line,1,last,cntr%anneri,erread) + 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 ELSEIF ( keyword_contains(previous_line,'ELECTRONS',alias='ELECTRON') ) THEN cntl%annee=.TRUE. CALL readsr(line,1,last,cntr%annere,erread) @@ -4034,6 +4049,36 @@ SUBROUTINE control locpot2%tlpot=.TRUE.! cmb-kk -Printing local potential ELSEIF ( keyword_contains(line,'DCACP') ) THEN vdwl%dcacp = .TRUE. + 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 ELSE ! Unknown keyword IF (' '/=line) THEN diff --git a/src/cotr.mod.F90 b/src/cotr.mod.F90 index 2ffcae7..3b11882 100644 --- a/src/cotr.mod.F90 +++ b/src/cotr.mod.F90 @@ -70,6 +70,9 @@ MODULE cotr INTEGER, ALLOCATABLE, SAVE :: lskptr(:,:) INTEGER, ALLOCATABLE, SAVE :: mm_ASKEL(:,:) + INTEGER, ALLOCATABLE, SAVE :: const_nconn(:) + INTEGER, ALLOCATABLE, SAVE :: const_conn(:,:) + REAL(real_8), ALLOCATABLE, SAVE :: hess(:,:) REAL(real_8), ALLOCATABLE, SAVE :: dtm(:) REAL(real_8), ALLOCATABLE, SAVE :: pmall(:) diff --git a/src/cpmd.F90 b/src/cpmd.F90 index d7f0235..63c6d33 100644 --- a/src/cpmd.F90 +++ b/src/cpmd.F90 @@ -97,6 +97,8 @@ SUBROUTINE cpmd USE cp_cuda_utils, ONLY: cp_cuda_init, cp_cuda_finalize USE ortho_utils, ONLY: ortho_init, ortho_finalize USE mts_utils, ONLY: read_mts_input + USE mimic_wrapper, ONLY: mimic_read_input, mimic_ifc_handshake,& + mimic_ifc_request_sizes, mimic_ifc_destroy IMPLICIT NONE CHARACTER(*), PARAMETER :: procedureN = 'cpmd' @@ -161,6 +163,13 @@ SUBROUTINE cpmd ! READ ATOMS, PSEUDOPOTENTIALS, COORDINATES CALL detsp ! detect number of species. + + IF (cntl%mimic) THEN + CALL mimic_read_input + CALL mimic_ifc_handshake + CALL mimic_ifc_request_sizes + endif + CALL mm_init ! set up some QM/MM stuff ! also needed for non-qmmm runs, so we can ! reduce use of #if defined(__GROMOS)/#endif @@ -410,7 +419,11 @@ SUBROUTINE cpmd call Delete(bicanonicalCpmdConfig) if (biCanonicalEnsembleDo) call Delete(bicanonicalCpmdInputConfig) - + + IF (cntl%mimic) THEN + CALL mimic_ifc_destroy + END IF + CALL finalize_cp_grp() CALL mp_end() ! ==--------------------------------------------------------------== diff --git a/src/detdof_utils.mod.F90 b/src/detdof_utils.mod.F90 index a9b7eb6..3ee11b4 100644 --- a/src/detdof_utils.mod.F90 +++ b/src/detdof_utils.mod.F90 @@ -91,7 +91,7 @@ SUBROUTINE detdof(tau0,tscr) glib=REAL(cotc0%nodim,kind=real_8) ELSE glib=REAL(cotc0%nodim,kind=real_8)-3._real_8 - IF (isos1%tisos.AND..NOT.lqmmm%qmmm) THEN + IF (isos1%tisos.AND.(.NOT.lqmmm%qmmm .AND..NOT.cntl%mimic)) THEN IF (ions1%nat.EQ.1) THEN glib=glib ELSEIF (ions1%nat.EQ.2) THEN @@ -167,7 +167,7 @@ SUBROUTINE detdof(tau0,tscr) ! == CONSTRAINTS == ! ==--------------------------------------------------------------== CALL dum2(tau0,tscr) - IF (cotc0%mcnstr.GT.0 .OR. cotr007%mrestr.GT.0) THEN + IF ((cotc0%mcnstr.GT.0 .OR. cotr007%mrestr.GT.0).AND..NOT.cntl%new_constraints) THEN l=MAX(cotc0%mcnstr,cotr007%mrestr) ALLOCATE(anorm(cotc0%nodim,l),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& @@ -177,61 +177,61 @@ SUBROUTINE detdof(tau0,tscr) __LINE__,__FILE__) ENDIF IF (cotc0%mcnstr.GT.0) 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__) ALLOCATE(xlagr(cotc0%nodim),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) ALLOCATE(ylagr(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) - CALL zeroing(xlagr)!,cotc0%mcnstr) CALL zeroing(ylagr)!,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 + 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 ENDIF IF (cotr007%mrestr.GT.0) THEN ALLOCATE(rskel(cotc0%nodim,cotr007%mrestr,12),STAT=ierr) @@ -329,7 +329,7 @@ SUBROUTINE detdof(tau0,tscr) IF (paral%io_parent)& WRITE(6,*)'THERE ARE QUITE A LOT OF FIXED ATOMS' ENDIF - IF (cotc0%mcnstr.GT.0.OR.cotr007%mrestr.GT.0) THEN + IF ((cotc0%mcnstr.GT.0.OR.cotr007%mrestr.GT.0).AND..NOT.cntl%new_constraints) THEN ALLOCATE(dxpar(3,maxsys%nax,maxsys%nsx),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) diff --git a/src/eextern_utils.mod.F90 b/src/eextern_utils.mod.F90 index 3adf268..b545488 100644 --- a/src/eextern_utils.mod.F90 +++ b/src/eextern_utils.mod.F90 @@ -11,6 +11,13 @@ MODULE eextern_utils USE kinds, ONLY: real_8 USE metr, ONLY: metr_com USE mm_input, ONLY: lqmmm + 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 USE parac, ONLY: parai,& paral USE pbc_utils, ONLY: pbc @@ -55,13 +62,36 @@ SUBROUTINE eextern(rhoe,v,eirop,tau0,fion,eext,tfor) CALL tiset(' EEXT',isub) ! ..ENERGY VEXT*RHOE eext=0._real_8 - IF (tfor.OR.cntl%texadd) THEN + 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 gfac=parm%omega/(spar%nr1s*spar%nr2s*spar%nr3s) eext=ddot(fpar%nnr1,rhoe,1,extf,1)*gfac + 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 ! ..ENERGY VEXT*RHOI AND GRADIENT eext2=0._real_8 eexti=0._real_8 - IF (.NOT.lqmmm%qmmm)THEN + IF (.NOT.lqmmm%qmmm.AND..NOT.cntl%mimic)THEN ! EGO-CPMD interface IF (cnti%iftype.EQ.1) THEN DO is=1,ions1%nsp diff --git a/src/finalp_utils.mod.F90 b/src/finalp_utils.mod.F90 index 7329a77..93d3bc7 100644 --- a/src/finalp_utils.mod.F90 +++ b/src/finalp_utils.mod.F90 @@ -74,7 +74,7 @@ SUBROUTINE finalp(tau0,fion,velp,eigv) ENDIF ENDIF CALL dumpr - CALL cnstpr + IF (.NOT.cntl%new_constraints) CALL cnstpr IF (cntl%tdiag) CALL wreigen(eigv,crge%f,ener_com%amu,crge%n) IF (paral%io_parent) WRITE(6,*) ! ==--------------------------------------------------------------== diff --git a/src/forcedr_driver.mod.F90 b/src/forcedr_driver.mod.F90 index 77e3ea6..290a7c8 100644 --- a/src/forcedr_driver.mod.F90 +++ b/src/forcedr_driver.mod.F90 @@ -140,7 +140,7 @@ SUBROUTINE forcedr(c0,c2,sc0,rhoe,psi,tau0,fion,eigv,& CALL restfc(tau0,fion) ENDIF !CALL mp_bcast(fion,3*maxsys%nax*maxsys%nsx,parai%source,parai%allgrp) - CALL mp_bcast(fion,3*maxsys%nax*maxsys%nsx,parai%io_source,parai%cp_grp) + CALL mp_bcast(fion,size(fion),parai%io_source,parai%cp_grp) ENDIF CALL mm_dim(mm_revert,oldstatus) CALL tihalt(procedureN,isub) diff --git a/src/forces_diag_utils.mod.F90 b/src/forces_diag_utils.mod.F90 index c122b39..21625b4 100644 --- a/src/forces_diag_utils.mod.F90 +++ b/src/forces_diag_utils.mod.F90 @@ -365,7 +365,7 @@ SUBROUTINE forces_diag(nstate,c0,c2,cr,sc0,cscr,vpp,eigv,& ! AK: NOTE: Restraints are already treated in FORCEDR() ENDIF !CALL mp_bcast(fion,3*maxsys%nax*maxsys%nsx,parai%source,parai%allgrp) - CALL mp_bcast(fion,3*maxsys%nax*maxsys%nsx,parai%io_source,parai%cp_grp) + CALL mp_bcast(fion,size(fion),parai%source,parai%cp_grp) ENDIF 200 CONTINUE ! Reset calste correctly even for soft exit diff --git a/src/forces_driver.mod.F90 b/src/forces_driver.mod.F90 index bceb1eb..de113eb 100644 --- a/src/forces_driver.mod.F90 +++ b/src/forces_driver.mod.F90 @@ -24,6 +24,7 @@ MODULE forces_driver USE jrotation_utils, ONLY: set_orbdist USE kinds, ONLY: real_8 USE kpts, ONLY: tkpts + USE mimic_wrapper, ONLY: mimic_energy USE mp_interface, ONLY: mp_sum USE nlforce_utils, ONLY: give_scr_nlforce,& nlforce @@ -234,6 +235,13 @@ SUBROUTINE forces(c0,c2,tau0,fion,rhoe,psi,nstate,nkpoint,lproj,tfor) ENDIF ener_com%exc=ener_com%exc+ehfx ener_com%etot=ener_com%etot+ehfx + 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 CALL tiset(procedureN//'_b',isub3) __NVTX_TIMER_START ( procedureN//'_b' ) @@ -422,7 +430,7 @@ SUBROUTINE forces(c0,c2,tau0,fion,rhoe,psi,nstate,nkpoint,lproj,tfor) 2000 CONTINUE rsactive = .FALSE. IF (tfor) THEN - CALL mp_sum(fion,3*maxsys%nax*maxsys%nsx,parai%allgrp) + CALL mp_sum(fion,size(fion),parai%allgrp) IF (paral%parent) THEN CALL symvec(fion) CALL taucl(fion) diff --git a/src/geofile_utils.mod.F90 b/src/geofile_utils.mod.F90 index f823f0e..2bac10c 100644 --- a/src/geofile_utils.mod.F90 +++ b/src/geofile_utils.mod.F90 @@ -13,6 +13,9 @@ MODULE geofile_utils USE linres, ONLY: tshl USE meta_multiple_walkers_utils, ONLY: mw_filename USE metr, ONLY: metr_com + USE mimic_wrapper, ONLY: mimic_revert_dim,& + mimic_save_dim,& + mimic_switch_dim USE mm_dim_utils, ONLY: mm_dim USE mm_dimmod, ONLY: clsaabox,& cpat,& @@ -63,6 +66,10 @@ SUBROUTINE geofile(tau0,velp,my_tag) ! ==--------------------------------------------------------------== tag=my_tag + IF (cntl%mimic) THEN + CALL mimic_save_dim() + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF CALL mm_dim(mm_go_mm,status) ALLOCATE(gr_tau(3,ions1%nat),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& @@ -239,6 +246,9 @@ SUBROUTINE geofile(tau0,velp,my_tag) WRITE(6,'(A,A)') ' GEOFILE| UNKNOWN TAG: ',tag ENDIF CALL mm_dim(mm_revert,status) + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF DEALLOCATE(gr_tau,STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'deallocation problem',& __LINE__,__FILE__) diff --git a/src/initrun_driver.mod.F90 b/src/initrun_driver.mod.F90 index 6ce8ff5..f2eb51a 100644 --- a/src/initrun_driver.mod.F90 +++ b/src/initrun_driver.mod.F90 @@ -13,6 +13,7 @@ MODULE initrun_driver taup,& velp USE copot_utils, ONLY: copot + USE cppt, ONLY: inyh USE elct, ONLY: crge USE error_handling, ONLY: stopgm USE fint, ONLY: fint1 @@ -23,6 +24,12 @@ MODULE initrun_driver USE linres, ONLY: lractive USE merge_utils, ONLY: merge USE metr, ONLY: metr_com + USE mimic_wrapper, ONLY: mimic_ifc_sort_fragments,& + mimic_revert_dim,& + mimic_save_dim,& + mimic_switch_dim,& + mimic_update_coords,& + mimic_control USE mm_dim_utils, ONLY: mm_dim USE mm_dimmod, ONLY: mm_go_mm,& mm_go_qm,& @@ -113,6 +120,9 @@ SUBROUTINE initrun(irec,c0,cm,sc0,rhoe,psi,eigv) CALL tiset(procedureN,isub) CALL mm_dim(mm_go_mm,status) + IF (cntl%mimic) THEN + CALL mimic_save_dim() + ENDIF IF (corel%tinlc) THEN ALLOCATE(vnlt(ncpw%nhg*2),STAT=ierr) @@ -146,14 +156,29 @@ SUBROUTINE initrun(irec,c0,cm,sc0,rhoe,psi,eigv) #endif IF (cntl%tshop) CALL stopgm('INITRUN','NEED A RESTART FILE',& __LINE__,__FILE__) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF ! Randomize cell parameters IF (cntl%tranc) CALL ranc(tau0) ! Randomization of the atomic coordinates. IF (cntl%tranp.AND.paral%io_parent) CALL ranp(tau0) - CALL mp_bcast(tau0,3*maxsys%nax*maxsys%nsx,parai%io_source,parai%cp_grp) + 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 ! INITIALIZATION OF WAVEFUNCTION CALL mm_dim(mm_go_qm,statusdummy) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF CALL rinitwf(c0,cm,sc0,nstate,tau0,taup,rhoe,psi) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ELSEIF (restart1%restart.AND.restart1%rnon) THEN IF (bsclcs.EQ.1) resthswf = .TRUE. #if defined (__GROMOS) @@ -161,21 +186,42 @@ SUBROUTINE initrun(irec,c0,cm,sc0,rhoe,psi,eigv) IF (rtr_l%restart_traj) CALL mm_restart_traj(tau0,velp,irec) ENDIF #endif + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF IF (restart1%rgeo.AND.paral%io_parent) CALL geofile(tau0,velp,'READ') + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF ! Randomize cell parameters IF (cntl%tranc) CALL ranc(tau0) ! Randomization of the atomic coordinates. IF (cntl%tranp.AND.paral%parent) CALL ranp(tau0) - CALL mp_bcast(tau0,3*maxsys%nax*maxsys%nsx,parai%io_source,parai%cp_grp) - CALL mp_bcast(velp,3*maxsys%nax*maxsys%nsx,parai%io_source,parai%cp_grp) + 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 ! INITIALIZATION OF WAVEFUNCTION CALL mm_dim(mm_go_qm,statusdummy) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF CALL rinitwf(c0,cm,sc0,nstate,tau0,taup,rhoe,psi) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF IF (restart1%rgeo.AND.paral%io_parent) CALL geofile(tau0,velp,'READ') - CALL mp_bcast(tau0,3*maxsys%nax*maxsys%nsx,parai%io_source,parai%cp_grp) - CALL mp_bcast(velp,3*maxsys%nax*maxsys%nsx,parai%io_source,parai%cp_grp) + CALL mp_bcast(tau0,size(tau0),parai%io_source,parai%cp_grp) + CALL mp_bcast(velp,size(velp),parai%io_source,parai%cp_grp) ELSEIF ((restart1%restart.AND. .NOT.restart1%rnon).OR.cntl%tmerge) THEN IF (bsclcs.EQ.1) resthswf = .FALSE. + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ! READ RESTART FILE IF (cntl%tmerge) THEN CALL MERGE(c0,cm,rhoe,psi) @@ -188,18 +234,31 @@ SUBROUTINE initrun(irec,c0,cm,sc0,rhoe,psi,eigv) IF (rtr_l%restart_traj) CALL mm_restart_traj(tau0,velp,irec) ENDIF #endif + 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 CALL mm_dim(mm_go_qm,statusdummy) CALL ainitwf(c0,nstate,tau0,taup,rhoe,psi) CALL mm_dim(mm_go_mm,statusdummy) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ! Randomize cell parameters IF (cntl%tranc) CALL ranc(tau0) ! Randomization of the atomic coordinates. IF (cntl%tranp.AND.paral%io_parent) CALL ranp(tau0) - CALL mp_bcast(tau0,3*maxsys%nax*maxsys%nsx,parai%io_source,parai%cp_grp) - CALL mp_bcast(velp,3*maxsys%nax*maxsys%nsx,parai%io_source,parai%cp_grp) + CALL mp_bcast(tau0,size(tau0),parai%io_source,parai%cp_grp) + CALL mp_bcast(velp,size(velp),parai%io_source,parai%cp_grp) ! Update symmetry operation IF (irec(irec_co).EQ.1) CALL updatsym(tau0) ! INITIALIZE A NEW WAVEFUNCTION + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF CALL mm_dim(mm_go_qm,statusdummy) IF (irec(irec_wf).EQ.0) THEN IF (cntl%tshop) CALL stopgm('INITRUN','NEED WAVEFUNCTIONS',& @@ -286,6 +345,9 @@ SUBROUTINE initrun(irec,c0,cm,sc0,rhoe,psi,eigv) trwanncx(ipx)=vdwwfl%trwannc ENDIF CALL mm_dim(mm_go_mm,statusdummy) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF IF (paral%io_parent) THEN IF (.NOT.lqmmm%qmmm) THEN ! we cannot do this since the gromos/cpmd translation @@ -317,6 +379,9 @@ SUBROUTINE initrun(irec,c0,cm,sc0,rhoe,psi,eigv) ! Initialization of density and potential ! for diagonalization schemes CALL mm_dim(mm_go_qm,statusdummy) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF IF (cntl%tdiag.AND..NOT.clc%classical) THEN IF (cntl%tshop) CALL stopgm('INITRUN','CNTL%TDIAG AND TSHOP NOT SUPPORTED'& ,& @@ -345,6 +410,9 @@ SUBROUTINE initrun(irec,c0,cm,sc0,rhoe,psi,eigv) iteropt%ndisrs=0 iteropt%ndistp=0 CALL mm_dim(mm_revert,status) + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF CALL tihalt(procedureN,isub) ! ==--------------------------------------------------------------== RETURN diff --git a/src/loadpa_utils.mod.F90 b/src/loadpa_utils.mod.F90 index 49ed615..d334a0c 100644 --- a/src/loadpa_utils.mod.F90 +++ b/src/loadpa_utils.mod.F90 @@ -16,6 +16,9 @@ MODULE loadpa_utils USE kinds, ONLY: real_8 USE kpts, ONLY: tkpts USE mm_dim_utils, ONLY: mm_dim + USE mimic_wrapper, ONLY: mimic_save_dim, & + mimic_switch_dim, & + mimic_revert_dim USE mm_dimmod, ONLY: mm_go_mm,& mm_revert USE mp_interface, ONLY: mp_sum,& @@ -27,7 +30,7 @@ MODULE loadpa_utils USE sphe, ONLY: gcutwmax USE system, ONLY: & fpar, iatpe, iatpt, ipept, mapgp, natpe, ncpw, nkpt, norbpe, parap, & - parm, spar + parm, spar, cntl USE timer, ONLY: tihalt,& tiset USE zeroing_utils, ONLY: zeroing @@ -104,6 +107,10 @@ SUBROUTINE loadpa ENDDO CALL mm_dim(mm_go_mm,oldstatus) + IF (cntl%mimic) THEN + CALL mimic_save_dim() + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ALLOCATE(iatpt(2,ions1%nat),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) @@ -115,6 +122,9 @@ SUBROUTINE loadpa iatpt(2,iat)=is ENDDO ENDDO + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF CALL mm_dim(mm_revert,oldstatus) DO i=0,parai%nproc-1 @@ -430,6 +440,11 @@ SUBROUTINE loadpa ! ==--------------------------------------------------------------== CALL leadim(parm%nr1,parm%nr2,parm%nr3,fpar%kr1,fpar%kr2,fpar%kr3) fpar%nnr1=fpar%kr1*fpar%kr2s*fpar%kr3s + + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF + DEALLOCATE(thread_buff,STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'deallocation problem', & __LINE__,__FILE__) diff --git a/src/md_driver.mod.F90 b/src/md_driver.mod.F90 index 6d43d4a..ad13a53 100644 --- a/src/md_driver.mod.F90 +++ b/src/md_driver.mod.F90 @@ -39,7 +39,10 @@ MODULE md_driver velp USE copot_utils, ONLY: copot,& give_scr_copot - USE cotr, ONLY: cotc0 + USE cotr, ONLY: cotc0,& + cnsval,& + ntcnst + USE cppt, ONLY: inyh USE ddipo_utils, ONLY: give_scr_ddipo USE detdof_utils, ONLY: detdof USE dispp_utils, ONLY: dispp @@ -99,6 +102,20 @@ MODULE md_driver USE meta_exlagr_utils, ONLY: ekincv_global USE meta_multiple_walkers_utils, ONLY: mw_assign_filenames,& mw_filename + 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 USE mfep, ONLY: mfepi USE mm_extrap, ONLY: cold , cold_high, & numcold, numcold_high, & @@ -168,7 +185,10 @@ MODULE md_driver tmpwr,& xf_fill_all_buffers,& xf_force_components - USE shake_utils, ONLY: cpmdshake + USE shake_utils, ONLY: cpmdshake,& + init_constraints,& + do_shake,& + do_rattle USE soc, ONLY: do_tsh_isc,& do_tsh_isc_lz,& isc_read_start,& @@ -249,6 +269,10 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) ! MTS] CALL tiset(procedureN,isub) + IF (cntl%mimic) THEN + CALL mimic_save_dim() + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF IF (cntl%tddft.AND.cntl%tresponse) CALL stopgm("MDDIAG",& "cntl%tddft.AND.cntl%tresponse NOT POSSIBLE",& __LINE__,__FILE__) @@ -288,6 +312,13 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) __LINE__,__FILE__) CALL zeroing(ch)!,nkpt%ngwk*n) ! CALL SETPOSOP + 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__) ELSE ALLOCATE(ch(1,1,1),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& @@ -380,12 +411,14 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) IF (fexist) CALL mp_bcast(ene_tri,SIZE(ene_tri),parai%io_source,parai%cp_grp) ENDIF ! SOC] - 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__) + 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 ! [EXACT FACTORIZATION IF (tshl%txfmqc) THEN ekincv=0._real_8 @@ -654,6 +687,10 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) IF (cnti%nomore.LT.0) GOTO 10000 ENDIF restf%nfnow=1 + IF (cntl%mimic) THEN + CALL mimic_ifc_init(rhoe, extf, tau0) + mimic_control%update_potential = .TRUE. + END IF ! ==--------------------------------------------------------------== ! TIME STEP FUNCTIONS CALL dynit(ekincp,ekin1,ekin2,temp1,temp2,ekinh1,ekinh2) @@ -671,6 +708,9 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) WRITE(6,*)& ' ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' ENDIF + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF ! ==--------------------------------------------------------------== ! == INITIALISATION == ! ==--------------------------------------------------------------== @@ -756,6 +796,9 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) CALL localize2(tau0,c0,c2,sc0,crge%n) ! ==--------------------------------------------------------------== ! + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF CALL mp_bcast(taup,SIZE(taup),parai%io_source,parai%cp_grp) IF (paral%parent) CALL dcopy(3*maxsys%nax*maxsys%nsx,taup,1,taui,1) ! ==--------------------------------------------------------------== @@ -768,6 +811,25 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) ! EHR] ! INITIALIZE VELOCITIES IF (paral%parent) CALL detdof(tau0,taur) + 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 ! INITIALIZE METADYNAMICS VARIABLES used also for ! SHOOTING from SADDLE POINT with RANDOM VELOCITIES @@ -781,11 +843,23 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) ener_com%erestr = 0.0_real_8 CALL rinvel(velp,c2,nstate) IF (paral%parent) CALL taucl(velp) - IF (paral%parent) CALL rattle(tau0,velp) + IF (paral%io_parent) THEN + IF (cntl%new_constraints) THEN + CALL do_rattle(tau0, velp) + ELSE + CALL rattle(tau0,velp) + END IF + END IF CALL rvscal(velp) ELSE IF (paral%parent) CALL taucl(velp) - IF (paral%parent) CALL rattle(tau0,velp) + IF (paral%io_parent) THEN + IF (cntl%new_constraints) THEN + CALL do_rattle(tau0, velp) + ELSE + CALL rattle(tau0,velp) + END IF + END IF IF (cntl%trescale) CALL rvscal(velp) ENDIF IF (cntl%quenchp) CALL zeroing(velp)!,3*maxsys%nax*maxsys%nsx) @@ -809,6 +883,9 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) +irec(irec_nop4) IF (cntl%tnosep .AND. itemp.EQ.0) CALL nospinit(ipwalk) ENDIF + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF ! CALL write_irec(irec) ! ..------------------------------------------------- @@ -872,6 +949,15 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) 'FORCES INITIALIZATION' WRITE(6,'(1X,64("="))') ENDIF + 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 ifcalc=0 IF (cntl%bsymm) THEN CALL bs_forces_diag(nstate,c0,c2,cm,sc0,cm(nx:),vpp,eigv,& @@ -929,6 +1015,10 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) IF (.NOT.cntl%bsymm) CALL dcopy(nnx,rin0,1,rm1(1),1) + IF (cntl%mimic) THEN + CALL mimic_sum_forces(fion) + END IF + IF (paral%io_parent) THEN CALL wrgeof(tau0,fion) CALL fileopen(3,filen,fo_app+fo_verb+fo_nofp,ferror) @@ -954,6 +1044,13 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) IF (paral%parent) CALL dispp(tau0,taui,disa) ! ENERGY OF THE NOSE THERMOSTATS CALL noseng(iteropt%nfi,velp,enose,enosp,dummy(1),1) + 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 econs=ekinp+ener_com%etot+enose+enosp+ener_com%ecnstr IF (cntl%cdft)econs=econs+cdftcom%cdft_v(1)*(cdftcom%cdft_nc+cdftcom%vgrad(1)) time2 =m_walltime() @@ -1007,6 +1104,9 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) ELSE hubbu%tpom=.False. ENDIF + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ! ANNEALING ! thermostats only if 1) standard dynamics 2) larger MTS step if ( .not.cntl%use_mts .or. mts_large_step ) then @@ -1033,7 +1133,11 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) ener_com%ecnstr = 0.0_real_8 IF (paral%io_parent) THEN CALL posupi(tau0,taup,velp) - IF (cotc0%mcnstr.NE.0) CALL cpmdshake(tau0,taup,velp) + IF (cntl%new_constraints) THEN + CALL do_shake(tau0, taup, velp) + ELSE + IF (cotc0%mcnstr.NE.0) CALL cpmdshake(tau0,taup,velp) + END IF #if defined (__QMECHCOUPL) IF (qmmech) THEN CALL mm_cpmd_update_links(taup, ions0%na, ions1%nsp, maxsys%nax, ions1%nat) @@ -1042,6 +1146,20 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) #endif ENDIF CALL mp_bcast(taup,SIZE(taup),parai%io_source,parai%cp_grp) + 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 ! Reset swap file to update (option BLOCK and ALL) IF (tkpts%tkblock) CALL reskpt_swap CALL phfac(taup) @@ -1077,7 +1195,7 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) ! switch off the info printing dmbi%inter_pt_firstcall=.FALSE. ELSE - + ifcalc = 0 IF (cntl%bsymm) THEN CALL bs_forces_diag(nstate,c0,c2,cm,sc0,cm(nx:),vpp,eigv,& rhoe,psi,& @@ -1175,6 +1293,9 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) ENDDO ENDIF ENDIF + IF (cntl%mimic) THEN + CALL mimic_sum_forces(fion) + END IF ! ================================================================== IF (ropt_mod%calste) CALL totstr @@ -1187,6 +1308,9 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) CALL mm_cpmd_elstat(taup,ions0%na,ions1%nsp,maxsys%nax,ions1%nat,c0,scr) ENDIF #endif + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ! FINAL UPDATE FOR VELOCITIES ener_com%ecnstr = 0.0_real_8 IF (paral%io_parent) THEN @@ -1197,9 +1321,13 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) ENDIF #endif IF (isos1%twall) CALL box_boundary(taup,velp) - CALL rattle(taup,velp) + IF (cntl%new_constraints) THEN + CALL do_rattle(taup, velp) + ELSE + CALL rattle(taup,velp) + END IF ENDIF - IF (paral%io_parent) CALL geofile(taup,velp,'WRITE') + IF (paral%parent.AND..NOT.cntl%mimic) CALL geofile(taup,velp,'WRITE') ! COMPUTE THE IONIC TEMPERATURE TEMPP IF (paral%io_parent) THEN CALL ekinpp(ekinp,velp) @@ -1209,6 +1337,9 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) ELSE tempp=ekinp*factem*2._real_8/glib ENDIF + IF (cntl%mimic) THEN + CALL mimic_subsystem_temperatures(velp) + END IF ! ..TSH[ tempp_sh=tempp ! ..TSH] @@ -1392,6 +1523,9 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) ENDIF ENDIF #endif + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ! temperature ramping CALL tempramp(temp1,temp2) ! UPDATE IONIC POSITIONS @@ -1402,6 +1536,9 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) CALL dcopy(nnx,rin0,1,rm1(1),1) CALL dcopy(nnx,rinp(1),1,rin0,1) ENDIF + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF ! STOP THE RUN IF THE USER HAS SET THE SIGNAL 30 ! ..TSH[ IF (tshl%tdtully) THEN @@ -1466,6 +1603,16 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) CALL hpm_stop('MD LOOP') #endif + 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 + IF (cntl%cdft)THEN IF (cntl%cdft_weight)CALL write_w(wdiff,"FINAL") ENDIF @@ -1687,6 +1834,9 @@ SUBROUTINE mddiag(c0,cm,c1,c2,sc0,vpp,gamx,gamy) __LINE__,__FILE__) ENDIF ! SOC] + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF CALL tihalt(procedureN,isub) ! ==--------------------------------------------------------------== END SUBROUTINE mddiag diff --git a/src/mdmain_utils.mod.F90 b/src/mdmain_utils.mod.F90 index 2f4600f..a926768 100644 --- a/src/mdmain_utils.mod.F90 +++ b/src/mdmain_utils.mod.F90 @@ -32,7 +32,10 @@ MODULE mdmain_utils velp USE copot_utils, ONLY: copot,& give_scr_copot - USE cotr, ONLY: cotc0 + USE cotr, ONLY: cnsval,& + cotc0,& + ntcnst + USE cppt, ONLY: inyh USE csize_utils, ONLY: csize USE ddipo_utils, ONLY: ddipo,& give_scr_ddipo @@ -41,6 +44,7 @@ MODULE mdmain_utils USE detdof_utils, ONLY: detdof USE dispp_utils, ONLY: dispp USE dynit_utils, ONLY: dynit + USE efld, ONLY: extf USE ekinpp_utils, ONLY: ekinpp USE elct, ONLY: crge USE ener, ONLY: chrg,& @@ -84,6 +88,20 @@ MODULE mdmain_utils USE meta_exlagr_utils, ONLY: ekincv_global USE meta_multiple_walkers_utils, ONLY: mw_assign_filenames,& mw_filename + 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 USE mfep, ONLY: mfepi USE mp_interface, ONLY: mp_bcast,& mp_sync @@ -140,7 +158,10 @@ MODULE mdmain_utils sample_wait USE setbsstate_utils, ONLY: setbsstate USE setirec_utils, ONLY: write_irec - USE shake_utils, ONLY: cpmdshake + USE shake_utils, ONLY: cpmdshake,& + init_constraints,& + do_shake,& + do_rattle USE soft, ONLY: soft_com USE spin, ONLY: spin_mod USE store_types, ONLY: & @@ -216,6 +237,10 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) REAL(real_8) :: mm_temp #endif CALL tiset(procedureN,isub) + IF (cntl%mimic) THEN + CALL mimic_save_dim() + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF time1 =m_walltime() ! Walker ID for multiple walker MTD. MTD part is in grandparent ipwalk=1 @@ -240,14 +265,18 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) ftmp=filenbs(i1:i2)//'.' CALL mw_filename(ftmp,filenbs,mfepi%istring) ENDIF - ! - 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__) + ! + 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 CALL zeroing(fion)!,SIZE(fion)) ALLOCATE(taui(3,maxsys%nax,maxsys%nsx),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& @@ -335,6 +364,10 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) IF (cnti%nomore.LT.0) GOTO 10000 ENDIF restf%nfnow=1 + IF (cntl%mimic) THEN + CALL mimic_ifc_init(rhoe, extf, tau0) + mimic_control%update_potential = .TRUE. + END IF ! ==--------------------------------------------------------------== ! TIME STEP FUNCTIONS CALL dynit(ekincp,ekin1,ekin2,temp1,temp2,ekinh1,ekinh2) @@ -344,6 +377,9 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) IF ((cntl%tnosee.OR.cntl%tnosep).AND.paral%parent) CALL nosepa(ipwalk,ipwalk) ! Dont symmetrize density cntl%tsymrho=.FALSE. + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF ! ==--------------------------------------------------------------== ! == INITIALIZATION == ! ==--------------------------------------------------------------== @@ -384,6 +420,10 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) IF (hfxc3%twscr.OR.(vdwl%vdwd.AND..NOT.vdwwfl%trwannc)) THEN CALL localize2(tau0,c0,c2,sc0,nstate) ENDIF + + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF ! ! NN: BROKEN SYMMETRY: QUENCHING TO BO SURFACE OF BS STATE IF (cntl%quenchb) THEN @@ -432,9 +472,34 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) ENDIF ENDIF + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF + ! INITIALIZE VELOCITIES IF (paral%parent) CALL detdof(tau0,taur) + ! 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 + ! INITIALIZE METADYNAMICS VARIABLES used also for ! SHOOTING from SADDLE POINT with RANDOM VELOCITIES IF (ionode .AND. (lmeta%lcolvardyn .OR. lmeta%lsadpnt)) THEN @@ -446,11 +511,23 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) ener_com%erestr = 0.0_real_8 CALL rinvel(velp,cm,crge%n) IF (paral%parent) CALL taucl(velp) - IF (paral%parent) CALL rattle(tau0,velp) + IF (paral%parent) THEN + IF (cntl%new_constraints) THEN + CALL do_rattle(tau0, velp) + ELSE + CALL rattle(tau0,velp) + END IF + END IF CALL rvscal(velp) ELSE IF (paral%parent) CALL taucl(velp) - IF (paral%parent) CALL rattle(tau0,velp) + IF (paral%parent) THEN + IF (cntl%new_constraints) THEN + CALL do_rattle(tau0, velp) + ELSE + CALL rattle(tau0,velp) + END IF + END IF IF (cntl%trescale) CALL rvscal(velp) ENDIF ! FOR NO RESTART THE VELOCITIES OF THE COEFFICIENTS OF THE @@ -477,6 +554,9 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) ! ..Initailize Molecular Mechanics subsystem: ! ..------------------------------------------------- ! .. + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF #if defined (__QMECHCOUPL) IF (paral%parent) THEN ! --------------------------------------------------------- @@ -526,6 +606,15 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) bsclcs=1 CALL setbsstate ENDIF + 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 CALL forcedr(c0(:,:,1),c2(:,:,1),sc0(:,:,1),rhoe,psi,& TAU0,FION,EIGV,NSTATE,1,.FALSE.,.TRUE.) ! STORE THE SPIN DENSITIES FOR PRINTING @@ -569,6 +658,10 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) lquench,lmetares,resetcv,ekinc,ekinp) ENDIF ENDIF + + IF (cntl%mimic) THEN + CALL mimic_sum_forces(fion) + END IF if (biCanonicalEnsembleDo) then call CpmdEnergiesGradients(bicanonicalCpmdConfig, ener_com%etot) @@ -681,6 +774,9 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) IF (.NOT.paral%parent) ropt_mod%prteig=.FALSE. ropt_mod%engpri=cprint%tprint.AND.MOD(iteropt%nfi-1,cprint%iprint_step).EQ.0 tstrng=cntl%tpmin.AND.MOD(iteropt%nfi,cnti%nomore).EQ.0 + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ! ANNEALING bsclcs=1 CALL berendsen(velp,cm(1,1,1),nstate,dummy,ekinc,0.0_real_8) @@ -726,7 +822,11 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) IF (paral%parent) THEN IF (cotc0%mcnstr.NE.0) THEN - CALL cpmdshake(tau0,taup,velp) + IF (cntl%new_constraints) THEN + CALL do_shake(tau0, taup, velp) + ELSE + CALL cpmdshake(tau0,taup,velp) + ENDIF ENDIF #if defined (__QMECHCOUPL) IF (qmmech) THEN @@ -737,6 +837,20 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) ENDIF CALL mp_bcast(taup,SIZE(taup),parai%io_source,parai%cp_grp) ! allgrp-cp_grp bugfix + 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 CALL phfac(taup) IF (corel%tinlc) CALL copot(rhoe,psi,ropt_mod%calste) ! BROKEN SYMMETRY WF UPDATED USING VELOCITY VERLET @@ -880,6 +994,9 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) ENDDO ENDIF + IF (cntl%mimic) THEN + CALL mimic_sum_forces(fion) + END IF IF (ropt_mod%calste) CALL totstr #if defined (__QMECHCOUPL) @@ -896,6 +1013,10 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) call CpmdEnergiesGradients(bicanonicalCpmdConfig, ener_com%etot) end if + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF + ! FINAL UPDATE FOR VELOCITIES ener_com%ecnstr = 0.0_real_8 ener_com%erestr = 0.0_real_8 @@ -907,7 +1028,13 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) ENDIF #endif IF (isos1%twall) CALL box_boundary(taup,velp) - CALL rattle(taup,velp) + 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 ENDIF ! UPDATE BROKEN_SYMMETRY_WAVEFUNCTION_VELOCITIES bsclcs=1 @@ -923,7 +1050,7 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) gamy(1,2),nstate) ENDIF - IF (paral%parent) CALL geofile(taup,velp,'WRITE') + IF (paral%parent.AND..NOT.cntl%mimic) CALL geofile(taup,velp,'WRITE') ! COMPUTE THE IONIC TEMPERATURE TEMPP IF (paral%parent) THEN CALL ekinpp(ekinp,velp) @@ -997,6 +1124,9 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) ELSE tempp=ekinp*factem*2._real_8/glib ENDIF + IF (cntl%mimic) then + CALL mimic_subsystem_temperatures(velp) + ENDIF ENDIF ! RESCALE ELECTRONIC VELOCITIES IF ((.NOT.cntl%bsymm).AND.cntl%tc)& @@ -1137,10 +1267,16 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) #endif ! Calculate EPR and EFG properties (not active for broken symmetry). IF (.NOT.cntl%bsymm) CALL epr_efg(rhoe,psi,iteropt%nfi) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ! temperature ramping CALL tempramp(temp1,temp2) ! UPDATE IONIC POSITIONS CALL dcopy(3*maxsys%nax*maxsys%nsx,taup(1,1,1),1,tau0(1,1,1),1) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF ! STOP THE RUN IF THE USER HAS SET THE SIGNAL 30 IF (soft_com%exsoft) GOTO 100 ! ================================================================== @@ -1148,6 +1284,15 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) ! ================================================================== ENDDO 100 CONTINUE + 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 IF (wannl%twann) THEN CALL ddipo(taup,c0(:,:,1),cm(:,:,1),c2(:,:,1),sc0,nstate,center) CALL forcedr(c0(:,:,1),c2(:,:,1),sc0(:,:,1),rhoe,psi,taup,fion,eigv,& @@ -1244,6 +1389,9 @@ SUBROUTINE mdmain(c0,cm,c2,sc0,gamx,gamy) ! NN: CLOSING THE FILE BS_ENG IF ((cntl%bsymm.AND.paral%parent).AND.paral%io_parent)& CALL fileclose(77) + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF ! ==--------------------------------------------------------------== CALL tihalt(procedureN,isub) END SUBROUTINE mdmain diff --git a/src/mimic_wrapper.mod.F90 b/src/mimic_wrapper.mod.F90 new file mode 100644 index 0000000..678e102 --- /dev/null +++ b/src/mimic_wrapper.mod.F90 @@ -0,0 +1,1557 @@ +!--------------------------------------------------------------------------------------------------- +! 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 diff --git a/src/mm_init_utils.mod.F90 b/src/mm_init_utils.mod.F90 index 4574f15..751f0ec 100644 --- a/src/mm_init_utils.mod.F90 +++ b/src/mm_init_utils.mod.F90 @@ -12,6 +12,10 @@ MODULE mm_init_utils USE ions, ONLY: ions0 USE isos, ONLY: isos1 USE kinds, ONLY: real_8 + USE mimic_wrapper, ONLY: mimic_save_dim,& + mimic_switch_dim,& + mimic_revert_dim,& + mimic_control USE mm_dimmod, ONLY: & clsaabox, cpat, cpsp, gratom, inc_l, mm_charge, mm_stat, mmdim, nam, & naq, nat_cpmd, nat_grm, ncag_l, solsolv, solvvv @@ -31,7 +35,8 @@ MODULE mm_init_utils paral USE rmas, ONLY: rmass USE store_types, ONLY: cprint - USE system, ONLY: maxsys + USE system, ONLY: cntl,& + maxsys USE zeroing_utils, ONLY: zeroing #include "sizeof.h" @@ -369,6 +374,10 @@ SUBROUTINE mm_compat_init INTEGER :: i, ia, ierr, is, NATm + IF (cntl%mimic) THEN + CALL mimic_save_dim() + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF mm_stat=.TRUE. ! we are (and stay) in QM dimensions. lqmmm%qmmm_verbose=.FALSE. ! no verbose QM/MM output. @@ -390,8 +399,13 @@ SUBROUTINE mm_compat_init natm=natm+ions0%na(is) ENDDO mmdim%natm = natm - mmdim%natq=mmdim%natm - mmdim%naxq=maxsys%nax + 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 IF (cprint%maxwriteatom.GT.mmdim%natm) cprint%maxwriteatom=mmdim%natm CALL mp_bcast(mmdim%natm,parai%io_source,parai%cp_grp) @@ -456,6 +470,9 @@ SUBROUTINE mm_compat_init CALL mp_bcast(cpsp,SIZE(cpsp),parai%io_source,parai%cp_grp) CALL mp_bcast(cpat,SIZE(cpat),parai%io_source,parai%cp_grp) CALL mp_bcast(gratom,SIZE(gratom),parai%io_source,parai%cp_grp) + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF RETURN END SUBROUTINE mm_compat_init diff --git a/src/mp_interface.mod.F90 b/src/mp_interface.mod.F90 index 73fac21..b15936e 100644 --- a/src/mp_interface.mod.F90 +++ b/src/mp_interface.mod.F90 @@ -13,6 +13,9 @@ MODULE mp_interface real_8 USE machine, ONLY: m_flush,& m_walltime +#ifdef __MIMIC + USE mimic_main, ONLY: mimic_init +#endif USE para_global, ONLY: para_buff_size,& para_stack_buff_size,& para_use_mpi_in_place @@ -303,6 +306,9 @@ SUBROUTINE mp_start __LINE__,__FILE__) ! an external interface has to set mp_comm_world mp_comm_world = MPI_COMM_WORLD +#ifdef __MIMIC + call mimic_init(mp_comm_world) +#endif ENDIF CALL mpi_errhandler_set ( mp_comm_world, MPI_ERRORS_RETURN, ierr ) IF (ierr/=0) CALL stopgm('MPI_INIT','mpi_errhandler_set',& diff --git a/src/mts_utils.mod.F90 b/src/mts_utils.mod.F90 index 289a86f..70ae0ab 100644 --- a/src/mts_utils.mod.F90 +++ b/src/mts_utils.mod.F90 @@ -94,6 +94,11 @@ subroutine read_mts_input write(output_unit,*) ' WARNING: MTS section will be ignored:', & ' USE_MTS keyword missing in CPMD section!' end if + ! 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 ! ! Main loop ! diff --git a/src/nosalloc_utils.mod.F90 b/src/nosalloc_utils.mod.F90 index f8a4cee..e1982dc 100644 --- a/src/nosalloc_utils.mod.F90 +++ b/src/nosalloc_utils.mod.F90 @@ -8,6 +8,9 @@ MODULE nosalloc_utils fo_verb USE ions, ONLY: ions0,& ions1 + USE mimic_wrapper, ONLY: mimic_save_dim,& + mimic_revert_dim,& + mimic_switch_dim USE mm_dim_utils, ONLY: mm_dim USE mm_dimmod, ONLY: mm_go_mm,& mm_revert @@ -58,6 +61,10 @@ SUBROUTINE nosalloc ! ==--------------------------------------------------------------== CALL tiset(procedureN,isub) + IF (cntl%mimic) THEN + CALL mimic_save_dim() + CALL mimic_switch_dim(go_qm=.FALSE.) + END IF ! initialize oldstatus to keep overly eager compilers to optimize it away. oldstatus=.FALSE. IF (lqmmm%qmmm)CALL mm_dim(mm_go_mm,oldstatus) @@ -286,6 +293,10 @@ SUBROUTINE nosalloc IF (lqmmm%qmmm)CALL mm_dim(mm_revert,oldstatus) + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + END IF + 120 FORMAT(i12,i12,f12.2,f12.2,f12.4 ) 121 FORMAT(i12,i12,i12,f12.2,f12.2) 122 FORMAT(a12,a12,a12,a12,a12) diff --git a/src/nosepa_utils.mod.F90 b/src/nosepa_utils.mod.F90 index 54eba5b..595f375 100644 --- a/src/nosepa_utils.mod.F90 +++ b/src/nosepa_utils.mod.F90 @@ -35,6 +35,7 @@ MODULE nosepa_utils USE system, ONLY: cnti,& cntl,& cntr,& + iatpt,& maxsp,& maxsys,& spar @@ -167,7 +168,7 @@ SUBROUTINE nosepa(np1,np2) IF (.NOT.(cntl%tpath.AND.cntl%tpimd)) THEN IF (ddof.LT.1.e-2_real_8) THEN dofst=dofst-3._real_8 - IF (isos1%tisos.AND..NOT.lqmmm%qmmm) THEN + IF (isos1%tisos.AND.(.NOT.lqmmm%qmmm.AND..NOT.cntl%mimic)) THEN IF (ions1%nat.EQ.1) THEN dofst=dofst ELSEIF (ions1%nat.EQ.2) THEN @@ -768,9 +769,13 @@ SUBROUTINE ilcttab 'RANGE START GREATER THAN RANGE STOP',& __LINE__,__FILE__) DO ia=istart,istop - ica=NAT_cpmd(ia) - isp=cpsp(ia) - ias=cpat(ia) + IF (cntl%mimic) THEN + ias=iatpt(1,ia) + isp=iatpt(2,ia) + ELSE + isp=cpsp(ia) + ias=cpat(ia) + END IF IF (lcttab(isp,ias).NE.0) THEN IF (lcttab(isp,ias).EQ.idxlct) THEN IF (paral%io_parent)& @@ -826,18 +831,16 @@ SUBROUTINE ilctmemb ENDDO ENDDO ! DEBUG - 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 - - + ! 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 ! ==--------------------------------------------------------------== RETURN END SUBROUTINE ilctmemb @@ -857,18 +860,31 @@ SUBROUTINE loccnst ! ==--------------------------------------------------------------== IF (cotc0%mcnstr.EQ.0) RETURN - ALLOCATE(n1(3*maxsys%nax*mmdim%nspm+1),STAT=ierr) + 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 IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) CALL zeroing(ddc)!,loct%nloct) iat=0 n1(1)=0 - DO is=1,mmdim%nspm - DO ia=1,NAm(is) - iat=iat+1 - n1(iat+1)=lcttab(is,ia) + 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 ENDDO - ENDDO + ELSE + DO is=1,mmdim%nspm + DO ia=1,NAm(is) + iat=iat+1 + n1(iat+1)=lcttab(is,ia) + ENDDO + ENDDO + END IF DO i=1,cotc0%mcnstr il1=n1(ntcnst(2,i)+1) ddc(il1)=ddc(il1)+1._real_8 @@ -895,7 +911,7 @@ SUBROUTINE loccnst WRITE(6,*) ' CONSTRAINTS ONLY ON EQUAL LOCAL THERMOSTATS ALLOWED ' IF (paral%io_parent)& WRITE(6,*) ' WITH THIS TYPE OF NOSE THERMOSTATS ' - CALL stopgm('NOSCNST','NOSE + CONSTRAINTS',& + CALL stopgm(procedureN,'NOSE + CONSTRAINTS',& __LINE__,__FILE__) END SUBROUTINE loccnst ! ================================================================== diff --git a/src/pbc_utils.mod.F90 b/src/pbc_utils.mod.F90 index 23144eb..fedbf50 100644 --- a/src/pbc_utils.mod.F90 +++ b/src/pbc_utils.mod.F90 @@ -5,6 +5,7 @@ MODULE pbc_utils USE isos, ONLY: isos1 USE kinds, ONLY: real_8 USE metr, ONLY: metr_com + USE mimic_wrapper, ONLY: mimic_control USE mm_dimmod, ONLY: clsaabox,& mm_stat USE system, ONLY: cntl, parm @@ -105,8 +106,14 @@ SUBROUTINE pbc(x1,y1,z1,x2,y2,z2,m,ap,ibrv) __LINE__,__FILE__) ENDIF ENDIF - ! ==--------------------------------------------------------------== - IF (isos1%tclust.AND.mm_stat) THEN + 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 + ! ==--------------------------------------------------------------== + IF (isos1%tclust.AND.mm_stat.AND..NOT.cntl%mimic) THEN IF (isos1%toned) THEN CALL pbc1(x1,y1,z1,x2,y2,z2,m,a,ibrv) ELSEIF (isos1%ttwod) THEN diff --git a/src/pcgrad_driver.mod.F90 b/src/pcgrad_driver.mod.F90 index 9c94acf..5da50a1 100644 --- a/src/pcgrad_driver.mod.F90 +++ b/src/pcgrad_driver.mod.F90 @@ -12,6 +12,10 @@ MODULE pcgrad_driver USE func, ONLY: func1 USE hubbardu, ONLY: hubbu USE kinds, ONLY: real_8 + USE mimic_wrapper, ONLY: mimic_save_dim,& + mimic_switch_dim,& + mimic_revert_dim,& + mimic_energy USE mm_dim_utils, ONLY: mm_dim USE mm_dimmod, ONLY: mm_go_mm,& mm_go_qm,& @@ -111,10 +115,17 @@ SUBROUTINE pcgrad(c0,c2,sc0,vpp,hnm1,rhoe,psi,tau0,nstate,dinit) ! we need to 'go mm' here so that the TSCR offset is correct. IF (lqmmm%qmmm) CALL mm_dim(mm_go_mm,oldstatus) ! ==--------------------------------------------------------------== + IF (cntl%mimic) THEN + CALL mimic_save_dim() + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ! TODO align for BG ALLOCATE(tscr(3,maxsys%nax,maxsys%nsx),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem', & __LINE__,__FILE__) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF ! ==--------------------------------------------------------------== IF (lqmmm%qmmm) CALL mm_dim(mm_go_qm,statusdummy) ! ==--------------------------------------------------------------== @@ -226,6 +237,9 @@ SUBROUTINE pcgrad(c0,c2,sc0,vpp,hnm1,rhoe,psi,tau0,nstate,dinit) ENDIF ENDIF IF (lqmmm%qmmm) CALL mm_dim(mm_revert,oldstatus) + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF DEALLOCATE(tscr,STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'deallocation problem', & __LINE__,__FILE__) @@ -541,6 +555,13 @@ SUBROUTINE xetot(nstate,x,sc0,tau0,tscr,rhoe,psi) CALL ortho(nstate,x,sc0) CALL rnlsm(x,nstate,1,1,.FALSE.) CALL rscpot(x,tau0,tscr,rhoe,psi,.FALSE.,.FALSE.,nstate,1) + 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 ENDIF ! ==--------------------------------------------------------------== RETURN diff --git a/src/potfor_utils.mod.F90 b/src/potfor_utils.mod.F90 index 95415ef..6c74496 100644 --- a/src/potfor_utils.mod.F90 +++ b/src/potfor_utils.mod.F90 @@ -50,11 +50,14 @@ SUBROUTINE potfor(fion,v,eirop) tyy, tzz, vcgs INTEGER :: ia, ig, ig1, is, isa, isub REAL(real_8) :: omtp + REAL(real_8), allocatable :: fion_temp(:,:,:) CALL tiset(procedureN,isub) omtp=2._real_8*parm%omega*parm%tpiba ig1=1 IF (geq0) ig1=2 + 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) #if defined (__VECTOR) IF (cntl%bigmem) THEN !$omp parallel do private(ISA,IA,IS,IG,EI123,RP,RHET,RHOG,RHETS,RHOGS, & @@ -77,9 +80,9 @@ SUBROUTINE potfor(fion,v,eirop) gz=CMPLX(0._real_8,gk(3,ig),kind=real_8) vcgs=scg(ig)*rhogs ei123=ei123*(rhops(is,ig)*vcgs+vps(is,ig)*rhets) - fion(1,ia,is)=fion(1,ia,is)+REAL(ei123*gx)*omtp - fion(2,ia,is)=fion(2,ia,is)+REAL(ei123*gy)*omtp - fion(3,ia,is)=fion(3,ia,is)+REAL(ei123*gz)*omtp + 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 ENDDO ENDDO ELSE @@ -104,9 +107,9 @@ SUBROUTINE potfor(fion,v,eirop) gz=CMPLX(0._real_8,gk(3,ig),kind=real_8) vcgs=scg(ig)*rhogs ei123=ei123*(rhops(is,ig)*vcgs+vps(is,ig)*rhets) - fion(1,ia,is)=fion(1,ia,is)+REAL(ei123*gx)*omtp - fion(2,ia,is)=fion(2,ia,is)+REAL(ei123*gy)*omtp - fion(3,ia,is)=fion(3,ia,is)+REAL(ei123*gz)*omtp + 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 ENDDO ENDDO ENDIF @@ -125,7 +128,7 @@ SUBROUTINE potfor(fion,v,eirop) !$omp parallel do private(IG,IS,IA,ISA) & !$omp private(RP,RHET,RHOG,RHETS,RHOGS,GX,GY,GZ,VCGS) & !$omp private(TXX,TYY,TZZ,EI123) & - !$omp reduction(+:FION) + !$omp reduction(+:FION_TEMP) #endif #endif DO ig=ig1,ncpw%nhg @@ -146,9 +149,9 @@ SUBROUTINE potfor(fion,v,eirop) DO ia=1,ions0%na(is) isa=isa+1 ei123=eigrb(ig,isa) - fion(1,ia,is)=fion(1,ia,is)+REAL(ei123*txx)*omtp - fion(2,ia,is)=fion(2,ia,is)+REAL(ei123*tyy)*omtp - fion(3,ia,is)=fion(3,ia,is)+REAL(ei123*tzz)*omtp + 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 ENDDO ENDDO ENDDO @@ -165,7 +168,7 @@ SUBROUTINE potfor(fion,v,eirop) !$omp parallel do private(IG,IS,IA,ISA) & !$omp private(RP,RHET,RHOG,RHETS,RHOGS,GX,GY,GZ,VCGS) & !$omp private(TXX,TYY,TZZ,EI123) & - !$omp reduction(+:FION) + !$omp reduction(+:FION_TEMP) #endif #endif DO ig=ig1,ncpw%nhg @@ -187,14 +190,15 @@ SUBROUTINE potfor(fion,v,eirop) isa=isa+1 ei123=ei1(isa,inyh(1,ig))*ei2(isa,inyh(2,ig))*& ei3(isa,inyh(3,ig)) - fion(1,ia,is)=fion(1,ia,is)+REAL(ei123*txx)*omtp - fion(2,ia,is)=fion(2,ia,is)+REAL(ei123*tyy)*omtp - fion(3,ia,is)=fion(3,ia,is)+REAL(ei123*tzz)*omtp + 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 ENDDO ENDDO ENDDO ENDIF #endif + fion(1:3, 1:maxval(ions0%na(1:ions1%nsp)), 1:ions1%nsp) = fion_temp(:,:,:) CALL tihalt(procedureN,isub) ! ==--------------------------------------------------------------== RETURN diff --git a/src/printp_utils.mod.F90 b/src/printp_utils.mod.F90 index c8fede5..b1f7455 100644 --- a/src/printp_utils.mod.F90 +++ b/src/printp_utils.mod.F90 @@ -177,7 +177,7 @@ SUBROUTINE printp(taur,taup,velp) CALL fileopen(4,f2,fo_app+if2,ferror) ENDIF IF (ropt_mod%rprint.OR.(ropt_mod%txyz.AND.rout1%xtout).OR.(ropt_mod%tdcd.AND.rout1%dcout)) THEN - IF (cotc0%mcnstr.GT.0.OR.cotr007%mrestr.GT.0) THEN + IF ((cotc0%mcnstr.GT.0.OR.cotr007%mrestr.GT.0).AND..NOT.cntl%new_constraints) THEN IF (paral%io_parent)& CALL fileopen(31,f3,if3,ferror) IF ((cotr007%mrestr.GT.0).AND.paral%io_parent)& @@ -280,7 +280,7 @@ SUBROUTINE printp(taur,taup,velp) IF (paral%io_parent)& CALL fileclose(32) ENDIF - IF (cotc0%mcnstr.GT.0.OR.cotr007%mrestr.GT.0) THEN + IF ((cotc0%mcnstr.GT.0.OR.cotr007%mrestr.GT.0).AND..NOT.cntl%new_constraints) THEN DO j=1,cotc0%mcnstr fval=fv(j) ityp=ntcnst(1,j) diff --git a/src/proppt_utils.mod.F90 b/src/proppt_utils.mod.F90 index faf9885..33b8dcc 100644 --- a/src/proppt_utils.mod.F90 +++ b/src/proppt_utils.mod.F90 @@ -23,6 +23,7 @@ MODULE proppt_utils rsdipo USE dipomod, ONLY: moment USE dist_prowfn_utils, ONLY: dist_prowfn + USE efld, ONLY: extf USE eicalc_utils, ONLY: eicalc USE elct, ONLY: crge USE elstpo_utils, ONLY: elstpo @@ -67,6 +68,16 @@ MODULE proppt_utils USE lodp, ONLY: & dmomlo, exd, extd, focc, nsdip, numbld, rcc, trmom, xmaxld, xminld, & ymaxld, yminld, zmaxld, zminld + 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 USE mm_dim_utils, ONLY: mm_dim USE mm_dimmod, ONLY: mm_go_mm,& mm_go_qm,& @@ -177,6 +188,10 @@ SUBROUTINE proppt IF ((cntl%tqmmm.AND.paral%parent).AND.paral%io_parent)& WRITE(6,*) 'PROPPT| WARNING QM/MM UNTESTED' + IF (cntl%mimic) THEN + CALL mimic_save_dim() + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ! ==--------------------------------------------------------------== IF (.NOT.soft_com%exsoft) CALL testex(soft_com%exsoft) IF (soft_com%exsoft) RETURN @@ -219,10 +234,15 @@ SUBROUTINE proppt IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) CALL mm_dim(mm_go_mm,oldstatus) - ALLOCATE(taup(3,maxsys%nax,maxsys%nsx),STAT=ierr) - IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& - __LINE__,__FILE__) + 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 CALL mm_dim(mm_go_qm,statusdummy) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF IF (cntl%tdiag.OR.coresl%tcores.OR.condpa%tconduct.OR.cldos%tldos) THEN ALLOCATE(rhoo(fpar%nnr1,clsd%nlsd),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& @@ -245,10 +265,16 @@ SUBROUTINE proppt ! ==--------------------------------------------------------------== IF (prop1%ceig) restart1%reigv=.TRUE. CALL read_irec(irec) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF CALL mm_dim(mm_go_mm,statusdummy) CALL zhrwf(1,irec,c0,cs,prop2%numorb,eigv,tau0,taup,taup,iteropt%nfi) CALL mp_bcast(tau0,SIZE(tau0),parai%io_source,parai%cp_grp) CALL mm_dim(mm_go_qm,statusdummy) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF ALLOCATE(fback(prop2%numorb),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) @@ -344,6 +370,9 @@ SUBROUTINE proppt ALLOCATE(eivps(ncpw%nhg),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) + IF (cntl%mimic) THEN + CALL mimic_ifc_init(rhoe, extf, tau0) + END IF ! ==------------------------------------------------------------== CALL phfac(tau0) CALL rnlsm(c0(:,:,1),prop2%numorb,1,1,.FALSE.) @@ -434,6 +463,9 @@ SUBROUTINE proppt ! == Calculates conductivity == ! ==--------------------------------------------------------------== IF (condpa%tconduct) THEN + IF (cntl%mimic) THEN + CALL mimic_ifc_init(rhoe, extf, tau0) + END IF CALL conductivity(c0,eigv,crge%f,crge%n,nkpt%nkpnt) IF (paral%io_parent) THEN ! cmb WRITE(6,'(A)')& @@ -459,6 +491,13 @@ SUBROUTINE proppt ALLOCATE(psi(il_psi_1d,il_psi_2d),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) + 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 CALL phfac(tau0) ! Compute electronic density in real space IF (tkpts%tkpnt) THEN @@ -482,6 +521,9 @@ SUBROUTINE proppt ! == Calculates core optical spectra == ! ==--------------------------------------------------------------== IF (coresl%tcores) THEN + IF (cntl%mimic) THEN + CALL mimic_ifc_init(rhoe, extf, tau0) + END IF CALL phfac(tau0) CALL core_spectra(tau0,c0,eigv,crge%f,crge%n,nkpt%nkpnt) ENDIF @@ -495,6 +537,9 @@ SUBROUTINE proppt ALLOCATE(sc0(nkpt%ngwk,nc/nkpt%ngwk),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) + IF (cntl%mimic) THEN + CALL mimic_ifc_init(rhoe, extf, tau0) + END IF CALL phfac(tau0) CALL localize(tau0,c0,cs,sc0,crge%n) CALL write_irec(irec) @@ -508,6 +553,9 @@ SUBROUTINE proppt ! ==--------------------------------------------------------------== IF (prop1%pwfn) THEN ! Orthogonalize + IF (cntl%mimic) THEN + CALL mimic_ifc_init(rhoe, extf, tau0) + END IF CALL phfac(tau0) IF (pslo_com%tivan) THEN CALL rnlsm(c0(:,:,1),prop2%numorb,1,1,.FALSE.) @@ -527,6 +575,9 @@ SUBROUTINE proppt ! ==--------------------------------------------------------------== IF (prop1%pylm) THEN lr=1 + IF (cntl%mimic) THEN + CALL mimic_ifc_init(rhoe, extf, tau0) + END IF CALL give_scr_ortho(lo,tag,prop2%numorb) lscr=MAX(lo,lr) CALL ortho(prop2%numorb,c0(:,:,1),cs) @@ -544,6 +595,9 @@ SUBROUTINE proppt ALLOCATE(psi(il_psi_1d,il_psi_2d),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) + IF (cntl%mimic) THEN + CALL mimic_ifc_init(rhoe, extf, tau0) + END IF CALL phfac(tau0) CALL ldos(c0,psi,eigv,crge%n) ENDIF @@ -572,6 +626,9 @@ SUBROUTINE proppt ALLOCATE(eivps(ncpw%nhg),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) + IF (cntl%mimic) THEN + CALL mimic_ifc_init(rhoe, extf, tau0) + END IF ! Compute electronic density in real space CALL phfac(tau0) IF (pslo_com%tivan) THEN @@ -679,6 +736,9 @@ SUBROUTINE proppt ALLOCATE(eivps(ncpw%nhg),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) + IF (cntl%mimic) THEN + CALL mimic_ifc_init(rhoe, extf, tau0) + END IF ! Compute electronic density in real space CALL phfac(tau0) IF (pslo_com%tivan) THEN @@ -801,6 +861,9 @@ SUBROUTINE proppt ! ==--------------------------------------------------------------== IF (prop1%glocl) THEN glocal%tgloc=.TRUE. + IF (cntl%mimic) THEN + CALL mimic_ifc_init(rhoe, extf, tau0) + END IF CALL phfac(tau0) IF (glocal%tg_real) THEN @@ -865,6 +928,9 @@ SUBROUTINE proppt ALLOCATE(eivps(nkpt%nhgk),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) + IF (cntl%mimic) THEN + CALL mimic_ifc_init(rhoe, extf, tau0) + END IF ! malloc complete. CALL phfac(tau0) CALL eicalc(eivps,eirop) @@ -1011,6 +1077,11 @@ SUBROUTINE proppt IF(ierr/=0) CALL stopgm(procedureN,'deallocation problem',& __LINE__,__FILE__) ENDIF + 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 ! ==--------------------------------------------------------------== DEALLOCATE(c0,STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'deallocation problem',& @@ -1024,6 +1095,9 @@ SUBROUTINE proppt CALL fnldealloc(.FALSE.,.FALSE.) ! ==--------------------------------------------------------------== CALL mm_dim(mm_revert,oldstatus) + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF RETURN END SUBROUTINE proppt ! ================================================================== diff --git a/src/ratom_utils.mod.F90 b/src/ratom_utils.mod.F90 index 271167c..187656a 100644 --- a/src/ratom_utils.mod.F90 +++ b/src/ratom_utils.mod.F90 @@ -36,6 +36,7 @@ MODULE ratom_utils USE meta_cell_utils, ONLY: meta_cell_inp USE meta_colvar_inp_utils, ONLY: meta_colvar_inp USE meta_dyn_def_utils, ONLY: meta_dyn_def + USE mimic_wrapper, ONLY: mimic_control USE mm_dimmod, ONLY: mmdim USE mm_input, ONLY: g96_vel,& lqmmm @@ -63,7 +64,8 @@ MODULE ratom_utils USE rmas, ONLY: rmass USE symm, ONLY: iun,& symmt - USE system, ONLY: maxsp,& + USE system, ONLY: cntl,& + maxsp,& maxsys USE tst2min_inp_utils, ONLY: tst2min_inp USE velocitinp_utils, ONLY: velocitinp @@ -159,6 +161,11 @@ SUBROUTINE ratom ! ==--------------------------------------------------------------== ! Initialization + IF (cntl%mimic) THEN + NSX_q = mimic_control%num_quantum_species + maxsys%nsx = mimic_control%num_species + maxsys%nax = mimic_control%max_atoms + ENDIF IF (lqmmm%qmmm)THEN NSX_q=mmdim%nspq maxsys%nsx=mmdim%nspm @@ -385,6 +392,7 @@ SUBROUTINE ratom ENDIF ! PSEUDOPOTENTIAL IF (lqmmm%qmmm)maxsys%nsx=mmdim%nspq + IF (cntl%mimic) maxsys%nsx = mimic_control%num_quantum_species IF (tnone) THEN CALL allelec(ions1%nsp,ecpnam) ELSEIF (tupf) THEN @@ -393,6 +401,7 @@ SUBROUTINE ratom CALL recpnew(ions1%nsp,ecpnam) ENDIF IF (lqmmm%qmmm)maxsys%nsx=mmdim%nspm + IF (cntl%mimic) maxsys%nsx = mimic_control%num_species IF (raggnew.GT.0._real_8) raggio(ions1%nsp)=raggnew ! ATOMIC COORDINATES IF (paral%io_parent)& @@ -748,6 +757,10 @@ SUBROUTINE ratom IF (lqmmm%qmmm)THEN maxsys%nsx=mmdim%nspq ENDIF + IF (cntl%mimic) THEN + maxsys%nsx = mimic_control%num_quantum_species + maxsys%nax = mimic_control%max_quantum_atoms + ENDIF ! ==--------------------------------------------------------------== RETURN END SUBROUTINE ratom diff --git a/src/rv30_utils.mod.F90 b/src/rv30_utils.mod.F90 index 1c0a3b0..3d87031 100644 --- a/src/rv30_utils.mod.F90 +++ b/src/rv30_utils.mod.F90 @@ -39,6 +39,9 @@ MODULE rv30_utils USE machine, ONLY: m_flush USE meta_multiple_walkers_utils, ONLY: mw_filename USE metr, ONLY: metr_com + USE mimic_wrapper, ONLY: mimic_save_dim,& + mimic_switch_dim,& + mimic_revert_dim USE mm_dimmod, ONLY: clsaabox USE mm_extrap, ONLY: cold,& nnow,& @@ -121,6 +124,10 @@ SUBROUTINE zhrwf(nr,irec,c0,cm,nstate,eigv,tau0,taum,taui,nfi) LOGICAL :: fexist CALL tiset(proceduren,isub) + IF (cntl%mimic) THEN + CALL mimic_save_dim() + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ! ==--------------------------------------------------------------== ! Construct the filename IF (paral%io_parent) THEN @@ -193,6 +200,9 @@ SUBROUTINE zhrwf(nr,irec,c0,cm,nstate,eigv,tau0,taum,taui,nfi) ' RESTART INFORMATION READ ON FILE ',filnam(ia:ie) CLOSE(nr) ENDIF + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF ! ==--------------------------------------------------------------== CALL tihalt(proceduren,isub) ! ==--------------------------------------------------------------== diff --git a/src/rwfopt_utils.mod.F90 b/src/rwfopt_utils.mod.F90 index 4fada5b..14cfd8a 100644 --- a/src/rwfopt_utils.mod.F90 +++ b/src/rwfopt_utils.mod.F90 @@ -16,10 +16,12 @@ MODULE rwfopt_utils give_scr_cplngs,& prcplngs USE cplngsmod + USE cppt, ONLY: inyh USE ddipo_utils, ONLY: give_scr_ddipo USE detdof_utils, ONLY: detdof USE dynit_utils, ONLY: dynit - USE efld + USE efld, ONLY: extf,& + textfld USE ehrenfest_utils, ONLY: ehrenfest USE elct USE enbandpri_utils, ONLY: enbandpri @@ -51,11 +53,23 @@ MODULE rwfopt_utils USE locpot USE lscal USE machine, ONLY: m_walltime + 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 USE mm_dim_utils, ONLY: mm_dim USE mm_dimmod USE mm_input USE mm_qmmm_forcedr_utils, ONLY: mm_qmmm_forcedr - USE mp_interface, ONLY: mp_sum + USE mp_interface, ONLY: mp_sum,& + mp_bcast USE norm USE parac, ONLY: parai,& paral @@ -78,7 +92,7 @@ MODULE rwfopt_utils USE syscomb_utils, ONLY: syscomb,& write_ksham USE system, ONLY: & - cnti, cntl, cntr, fpar, locpot2, maxsys, nacc, ncpw, nkpt, parm, spar + cnti, cntl, cntr, fpar, locpot2, maxsys, nacc, ncpw, nkpt, parm, spar, parap USE td_input USE td_utils, ONLY: getnorm_k,& load_ex_states,& @@ -148,9 +162,13 @@ SUBROUTINE rwfopt(c0,c2,sc0,pme,gde,vpp,eigv) tscr(:,:,:), vpotx3s(:,:,:), & wdsave(:) REAL(real_8), POINTER :: vpotx3(:,:,:) + COMPLEX(real_8), DIMENSION(:,:,:), ALLOCATABLE :: dummy CALL tiset(procedureN,isub) time1 =m_walltime() + IF (cntl%mimic) THEN + CALL mimic_save_dim() + ENDIF ! ==--------------------------------------------------------------== ! SCR creation CALL rhoe_psi_size(il_psi_1d=il_psi_1d,il_psi_2d=il_psi_2d,& @@ -190,9 +208,24 @@ SUBROUTINE rwfopt(c0,c2,sc0,pme,gde,vpp,eigv) ENDIF ENDIF + 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 + CALL mm_dim(mm_go_mm,oldstatus) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ! CB: TAUP,FION,VELP handled by BS_WFO - IF (.NOT.cntl%bsymm) THEN + IF (.NOT.cntl%bsymm.AND..NOT.cntl%mimic) THEN ALLOCATE(taup(3,maxsys%nax,maxsys%nsx),STAT=ierr) IF(ierr/=0) CALL stopgm(procedureN,'allocation problem',& __LINE__,__FILE__) @@ -248,6 +281,9 @@ SUBROUTINE rwfopt(c0,c2,sc0,pme,gde,vpp,eigv) ENDIF ! CALL mm_dim(mm_go_qm,statusdummy) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF nacc = 7 IF (cntl%tdiag.OR.cntl%tdiagopt) THEN nnx=fpar%nnr1*clsd%nlsd @@ -310,6 +346,15 @@ SUBROUTINE rwfopt(c0,c2,sc0,pme,gde,vpp,eigv) CALL write_ksham(c0,c2,sc0,rhoe,psi,eigv) GOTO 150 ENDIF + 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 ! McB this is not necessary in standard runs ! McB (maybe for cntl%cdft, but it is WRONG! for QM/MM (cf. setsys.F)!) ! IF(TCLUST.AND.TCENT.AND..NOT.TCLAS)CALL QM_CENTRE @@ -317,7 +362,13 @@ SUBROUTINE rwfopt(c0,c2,sc0,pme,gde,vpp,eigv) IF (paral%parent) THEN ! McB DETDOF needs to be called in full system dimensions CALL mm_dim(mm_go_mm,statusdummy) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF CALL detdof(tau0,tscr) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF CALL mm_dim(mm_go_qm,statusdummy) ENDIF IF (cntl%cdft)THEN @@ -499,6 +550,9 @@ SUBROUTINE rwfopt(c0,c2,sc0,pme,gde,vpp,eigv) update_pot=.FALSE. ENDIF CALL mm_dim(mm_go_qm,statusdummy) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF IF (cntl%tdiag) THEN CALL updrho(c0,c2,gde,sc0,pme,vpp,tau0,fion,eigv,& rhoe,psi,& @@ -521,6 +575,9 @@ SUBROUTINE rwfopt(c0,c2,sc0,pme,gde,vpp,eigv) ENDIF ENDIF CALL mm_dim(mm_go_mm,statusdummy) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF IF (paral%io_parent) THEN ropt_mod%engpri=MOD(infi-1,cprint%iprint_step).EQ.0 ELSE @@ -660,9 +717,15 @@ SUBROUTINE rwfopt(c0,c2,sc0,pme,gde,vpp,eigv) ENDIF ENDIF CALL mm_dim(mm_go_qm,statusdummy) + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF IF (rout1%rhoout.AND.(.NOT.cntl%bsymm.OR.ropt_mod%convwf)) THEN CALL rhopri(c0,tau0,rhoe,psi(:,1),crge%n,nkpt%nkpnt) ENDIF + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ! EHR[ ! ==--------------------------------------------------------------== ! == COMPUTE SPECTRA USING PROPAGATION OF PERTURBED WFs @@ -828,6 +891,9 @@ SUBROUTINE rwfopt(c0,c2,sc0,pme,gde,vpp,eigv) ENDIF ! Calculate ionic forces IF (ropt_mod%convwf.AND.(tfor.OR.cntl%tpres.OR.vdwl%vdwd)) THEN + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF IF (vdwl%vdwd) CALL localize(tau0,c0,c2,sc0,crge%n) IF (cntl%tdiag) THEN CALL updrho(c0,c2,gde,sc0,pme,vpp,tau0,fion,eigv,& @@ -845,12 +911,23 @@ SUBROUTINE rwfopt(c0,c2,sc0,pme,gde,vpp,eigv) rhoe,psi,crge%n,.TRUE.,update_pot) ENDIF ENDIF + IF (cntl%mimic) THEN + CALL mimic_sum_forces(fion) + ENDIF ! Calculate norm of nuclear gradient CALL gsize(fion,gnmax,gnorm) IF (cntl%tpres) THEN CALL totstr ENDIF + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ELSE + 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 gnmax=0.0_real_8 gnorm=0.0_real_8 ENDIF @@ -859,6 +936,9 @@ SUBROUTINE rwfopt(c0,c2,sc0,pme,gde,vpp,eigv) IF (lspin2%teprof.AND.lspin2%tlse) THEN CALL lseprof(c0,rhoe,psi,tau0,fion,crge%n) ENDIF + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ! EHR[ IF (cntl%tpspec) THEN IF (rout1%rhoout) CALL rhopri(c0,tau0,rhoe,psi(:,1),crge%n,nkpt%nkpnt) @@ -971,6 +1051,9 @@ SUBROUTINE rwfopt(c0,c2,sc0,pme,gde,vpp,eigv) IF(cntl%thubb)DEALLOCATE(c2u) ! ==--------------------------------------------------------------== CALL mm_dim(mm_revert,oldstatus) + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF CALL tihalt(procedureN,isub) RETURN END SUBROUTINE rwfopt diff --git a/src/setirec_utils.mod.F90 b/src/setirec_utils.mod.F90 index 9ea7e6c..7a746fd 100644 --- a/src/setirec_utils.mod.F90 +++ b/src/setirec_utils.mod.F90 @@ -121,7 +121,7 @@ SUBROUTINE read_irec(irec) irec(irec_prfo) = isetone(restart1%rlscl.AND.cntl%prfo) irec(irec_rdtl) = isetone(restart1%rdtl) irec(irec_cons) = isetone(restart1%rcon) - irec(irec_ctrans) = isetone(cntl%tqmmm) + irec(irec_ctrans) = isetone(cntl%tqmmm.OR.cntl%mimic) irec(irec_xtrp) = isetone(cntl%textrap.AND.restart1%rxtp) irec(irec_shop) = isetone(sh03%tshopres.AND.(.NOT.cntl%tfusi)) irec(irec_shopbo) = isetone(sh03%tshopres.AND.cntl%tmdbo.AND.(.NOT.cntl%tfusi)) @@ -190,7 +190,7 @@ SUBROUTINE write_irec(irec) irec(irec_prfo) = isetone(cntl%prfo) irec(irec_rdtl) = isetone(cntr%tolene.NE.0.0_real_8.OR.cntr%tolad.NE.0.0_real_8) irec(irec_cons) = isetone(cotc0%mcnstr.NE.0) - irec(irec_ctrans) = isetone(cntl%tqmmm) + irec(irec_ctrans) = isetone(cntl%tqmmm.OR.cntl%mimic) irec(irec_xtrp) = isetone(cntl%textrap.AND.cntl%tstrxtp) irec(irec_shop) = isetone(cntl%tshop.AND.(.NOT.cntl%tfusi)) irec(irec_shopbo) = isetone(cntl%tshop.AND.cntl%tmdbo.AND.(.NOT.cntl%tfusi)) diff --git a/src/setsys_utils.mod.F90 b/src/setsys_utils.mod.F90 index e5974f8..286f0f2 100644 --- a/src/setsys_utils.mod.F90 +++ b/src/setsys_utils.mod.F90 @@ -54,6 +54,9 @@ MODULE setsys_utils USE kpts, ONLY: kpts_com,& tkpts USE metr, ONLY: metr_com + USE mimic_wrapper, ONLY: mimic_save_dim,& + mimic_switch_dim,& + mimic_revert_dim USE mm_dim_utils, ONLY: mm_dim USE mm_dimmod, ONLY: clsaabox,& mm_go_mm,& @@ -159,6 +162,10 @@ SUBROUTINE setsys mmdim%nspm =maxsys%nsx ENDIF IF (lqmmm%qmmm)CALL mm_dim(mm_go_mm,status) + IF (cntl%mimic) THEN + CALL mimic_save_dim() + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF IF (.NOT.paral%io_parent) GOTO 9999 IF (cntl%tscale) THEN ! ==--------------------------------------------------------------== @@ -355,7 +362,8 @@ SUBROUTINE setsys clsaabox%mm_c_trans(3)=0.0_real_8 ! EHR[ IF (isos1%tclust.AND.isos1%tcent.AND..NOT.tclas& - .AND..NOT.cntl%tmdeh.AND..NOT.cntl%tpspec) THEN + .AND..NOT.cntl%tmdeh.AND..NOT.cntl%tpspec& + .AND..NOT.cntl%mimic) THEN ! EHR] IF (lqmmm%qmmm)THEN #if defined (__GROMOS) @@ -425,6 +433,9 @@ SUBROUTINE setsys IF (lqmmm%qmmm)THEN IF ( mmdim%natm.GT.500) GOTO 111 ENDIF + IF (cntl%mimic) THEN + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF WRITE(6,'(/,1X,29("*"),A,28("*"))') ' ATOMS ' WRITE(6,'(A,A)') ' NR TYPE X(BOHR) Y(BOHR)',& ' Z(BOHR) MBL' @@ -1620,6 +1631,9 @@ SUBROUTINE setsys CALL mp_bcast_byte(clsaabox, size_in_bytes_of(clsaabox), parai%io_source, parai%cp_grp) ! ==--------------------------------------------------------------== IF (lqmmm%qmmm)CALL mm_dim(mm_go_qm,status) ! because it is needed + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF CALL cpl_para RETURN END SUBROUTINE setsys diff --git a/src/shake_utils.mod.F90 b/src/shake_utils.mod.F90 index 25bf15a..3062554 100644 --- a/src/shake_utils.mod.F90 +++ b/src/shake_utils.mod.F90 @@ -3,6 +3,8 @@ MODULE shake_utils USE cotr, ONLY: & anorm, cnsval, cnsval_dest, cotc0, dtm, fc, fv, grate, kshove, & mm_askel, ntcnst, xlagr, ylagr + USE constraint + USE constraint_utils USE error_handling, ONLY: stopgm USE ions, ONLY: ions0,& ions1 @@ -23,12 +25,62 @@ MODULE shake_utils PRIVATE + type(constraint_matrix) :: const_matrix + PUBLIC :: cpmdshake + PUBLIC :: init_constraints + PUBLIC :: do_shake + PUBLIC :: do_rattle !public :: solvs !public :: andforces CONTAINS + 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 + ! ================================================================== SUBROUTINE cpmdshake(tau0,taup,velp) ! ==--------------------------------------------------------------== diff --git a/src/sysin_utils.mod.F90 b/src/sysin_utils.mod.F90 index 43e822c..bb11d9f 100644 --- a/src/sysin_utils.mod.F90 +++ b/src/sysin_utils.mod.F90 @@ -1445,6 +1445,13 @@ SUBROUTINE check_options() ! ENDIF + 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 + ! ! AK: sanity check. ibrav=-2 is used for implicit triclinic symmetry ! with given cell vectors. symmetry should not be used in that case. diff --git a/src/system.mod.F90 b/src/system.mod.F90 index 4eb8d3f..6e527c1 100644 --- a/src/system.mod.F90 +++ b/src/system.mod.F90 @@ -563,6 +563,10 @@ MODULE system LOGICAL :: thubb = .FALSE. LOGICAL :: use_mts = .FALSE. LOGICAL :: use_scaled_hfx = .FALSE. + LOGICAL :: mimic = .FALSE. + LOGICAL :: new_constraints = .FALSE. + LOGICAL :: pbicgstab = .FALSE. + LOGICAL :: anneal_dual = .FALSE. END TYPE cntl_t TYPE(cntl_t), SAVE, PUBLIC :: cntl ! ================================================================== @@ -673,6 +677,8 @@ MODULE system INTEGER :: isocs = HUGE(0) !vw not initialized at all INTEGER :: jsoct = HUGE(0) !vw not initialized at all INTEGER :: disortho_bsize = HUGE(0) + INTEGER :: shake_maxstep = HUGE(0) + INTEGER :: shake_cg_iter = HUGE(0) END TYPE cnti_t TYPE(cnti_t), SAVE, PUBLIC :: cnti ! ================================================================== @@ -790,6 +796,7 @@ MODULE system REAL(real_8) :: dampge = HUGE(0.0_real_8) REAL(real_8) :: dampgc = HUGE(0.0_real_8) REAL(real_8) :: gfreq = HUGE(0.0_real_8) !vw not initialized at all + REAL(real_8), DIMENSION(2) :: anneal_factors = HUGE(0.0_real_8) END TYPE cntr_t TYPE(cntr_t), SAVE, PUBLIC :: cntr ! strings diff --git a/src/updrho_utils.mod.F90 b/src/updrho_utils.mod.F90 index a5f79c4..de119a1 100644 --- a/src/updrho_utils.mod.F90 +++ b/src/updrho_utils.mod.F90 @@ -53,6 +53,7 @@ MODULE updrho_utils USE kpnt, ONLY: wk USE kpts, ONLY: tkpts USE ksdiag_utils, ONLY: ksdiag + USE mimic_wrapper, ONLY: mimic_energy USE mixing_g_utils, ONLY: give_scr_mixing,& mixing_g USE mixing_r_utils, ONLY: mixing_r @@ -655,7 +656,14 @@ SUBROUTINE updrho(c0,c2,cr,sc0,cscr,vpp,tau0,fion,eigv,rhoe,& ! CALL RNLRH(ENL,NSTATE,NKPTS) ener_com%etot=ener_com%eeig+(ener_com%exc-ener_com%vxc)+ener_com%eht& +vdwr%evdw ! Empirical van der Waals correction - ! Check total energy difference + 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 detot=ener_com%etot-etot0 IF (ABS(detot).LT.cntr%toldetot) THEN netot=netot+1 diff --git a/src/wrener_utils.mod.F90 b/src/wrener_utils.mod.F90 index cfb6b32..eca7f70 100644 --- a/src/wrener_utils.mod.F90 +++ b/src/wrener_utils.mod.F90 @@ -25,6 +25,8 @@ MODULE wrener_utils USE kpts, ONLY: tkpts USE machine, ONLY: m_flush USE metr, ONLY: metr_com + USE mimic_wrapper, ONLY: mimic_control,& + mimic_energy USE mm_dim_utils, ONLY: mm_dim USE mm_dimmod, ONLY: mm_go_qm,& mm_revert @@ -329,7 +331,9 @@ SUBROUTINE wrprint_md(eigv,f,amu,nstate,tau0,fion,& CALL wrgeo(tau0) ENDIF IF (ok_print(cprint%iprint(iprint_info),engpri,tend)) THEN - CALL cnstpr + IF (.NOT.cntl%new_constraints) THEN + CALL cnstpr + END IF IF (.NOT.cntl%bsymm) CALL wrener IF (.NOT.(infi.EQ.cnti%nomore.OR.convergence)) THEN IF (cntl%tmdbo) THEN @@ -539,6 +543,7 @@ SUBROUTINE wrener CHARACTER(len=14) :: defenergy, defenergy1, & defenergy2 CHARACTER(len=16) :: defenergy3 + CHARACTER(len=25) :: defenergy_mimic INTEGER :: is LOGICAL :: statusdummy REAL(real_8) :: cdft_ener, elec1, elec2, & @@ -683,6 +688,18 @@ SUBROUTINE wrener WRITE(6,'(1X,A,T22,A,T41,F20.8,A)')& ' ','DIFFERENCE =',& ener_com%etot-(eqm_r%eqm+eqm_r%eqmmm+eqm_r%EMM),' A.U.' + 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 ELSE IF (paral%io_parent)& WRITE(6,'(/,1X,A,T25,A,T41,F20.8,A)')& diff --git a/src/wrgeo_utils.mod.F90 b/src/wrgeo_utils.mod.F90 index 49e1483..a7b0897 100644 --- a/src/wrgeo_utils.mod.F90 +++ b/src/wrgeo_utils.mod.F90 @@ -9,6 +9,9 @@ MODULE wrgeo_utils USE ions, ONLY: ions0,& ions1 USE kinds, ONLY: real_8 + USE mimic_wrapper, ONLY: mimic_save_dim,& + mimic_switch_dim,& + mimic_revert_dim USE mm_dim_utils, ONLY: mm_dim USE mm_dimmod, ONLY: mm_go_mm,& mm_revert,& @@ -21,7 +24,8 @@ MODULE wrgeo_utils USE store_types, ONLY: cprint,& iprint_force,& rout1 - USE system, ONLY: cnti + USE system, ONLY: cnti,& + cntl USE timer, ONLY: tihalt,& tiset @@ -43,6 +47,10 @@ SUBROUTINE wrgeo(tau0) INTEGER :: ia, iat, is, k + IF (cntl%mimic) THEN + CALL mimic_save_dim() + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF IF (paral%io_parent)& WRITE(6,'(/,1X,64("*"))') IF (paral%io_parent)& @@ -74,6 +82,9 @@ SUBROUTINE wrgeo(tau0) IF (paral%io_parent)& WRITE(6,'(1X,64("*"),/)') ! ==--------------------------------------------------------------== + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF RETURN END SUBROUTINE wrgeo ! ================================================================== @@ -87,8 +98,16 @@ SUBROUTINE wrgeof(tau0,fion) CALL tiset(procedureN,isub) + IF (cntl%mimic) THEN + CALL mimic_save_dim() + CALL mimic_switch_dim(go_qm=.TRUE.) + ENDIF + IF (cprint%iprint(iprint_force).LE.0) THEN CALL wrgeo(tau0) + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF CALL tihalt(procedureN,isub) RETURN ENDIF @@ -120,6 +139,9 @@ SUBROUTINE wrgeof(tau0,fion) ENDDO ENDIF ! ==--------------------------------------------------------------== + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF CALL tihalt(procedureN,isub) RETURN END SUBROUTINE wrgeof @@ -161,6 +183,10 @@ SUBROUTINE wrgeox(tau0) LOGICAL :: ferror, status REAL(real_8), ALLOCATABLE :: gr_tau(:,:) + IF (cntl%mimic) THEN + CALL mimic_save_dim() + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF IF (rout1%xgout .AND. MOD(infi,cnti%ngxyz).EQ.0) THEN iunit=1 CALL mm_dim(mm_go_mm,status) @@ -204,6 +230,9 @@ SUBROUTINE wrgeox(tau0) IF(ierr/=0) CALL stopgm(procedureN,'deallocation problem',& __LINE__,__FILE__) ENDIF + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF ! ==--------------------------------------------------------------== RETURN END SUBROUTINE wrgeox diff --git a/src/wv30_utils.mod.F90 b/src/wv30_utils.mod.F90 index d1f9c5e..7eda6a8 100644 --- a/src/wv30_utils.mod.F90 +++ b/src/wv30_utils.mod.F90 @@ -42,6 +42,9 @@ MODULE wv30_utils USE machine, ONLY: m_system USE meta_multiple_walkers_utils, ONLY: mw_filename USE metr, ONLY: metr_com + USE mimic_wrapper, ONLY: mimic_save_dim,& + mimic_switch_dim,& + mimic_revert_dim USE mm_dimmod, ONLY: clsaabox USE mm_extrap, ONLY: cold,& nnow,& @@ -125,6 +128,10 @@ SUBROUTINE zhwwf (nw,irec,c0,cm,nstate,eigv,tau0,taum,taui,nfi) IF (ibench(1).EQ.1) RETURN ! ==--------------------------------------------------------------== CALL tiset(procedureN,isub) + IF (cntl%mimic) THEN + CALL mimic_save_dim() + CALL mimic_switch_dim(go_qm=.FALSE.) + ENDIF ! ==--------------------------------------------------------------== ! Construct the filename @@ -226,6 +233,9 @@ SUBROUTINE zhwwf (nw,irec,c0,cm,nstate,eigv,tau0,taum,taui,nfi) IF (paral%parent)THEN IF (cntl%cdft)CALL wcdft_restart() ENDIF + IF (cntl%mimic) THEN + CALL mimic_revert_dim() + ENDIF ! ==--------------------------------------------------------------== CALL tihalt(procedureN,isub) ! ==--------------------------------------------------------------==