!*********************************************************************** ! ! Copyright (c) 2015, Lawrence Livermore National Security, LLC. ! Produced at the Lawrence Livermore National ! Laboratory. ! Written by Nicolas Schunck, schunck1@llnl.gov ! ! LLNL-CODE-470611 All rights reserved. ! ! Copyright 2011, N. Schunck, J. Dobaczewski, J. McDonnell, ! W. Satula, J.A. Sheikh, A. Staszczak,A ! M. Stoitsov, P. Toivanen ! Copyright 2009, J. Dobaczewski, W. Satula, B.G. Carlsson, J. Engel, ! P. Olbratowski, P. Powalowski, M. Sadziak, ! J. Sarich, N. Schunck, A. Staszczak, M. Stoitsov, ! M. Zalewski, H. Zdunczuk ! Copyright 2004, 2005, J. Dobaczewski, P. Olbratowski ! Copyright 1997, 2000, J. Dobaczewski, J. Dudek ! ! This file is part of HFODD. ! ! HFODD is free software: you can redistribute it and/or modify it ! under the terms of the GNU General Public License as published by ! the Free Software Foundation, either version 3 of the License, or ! (at your option) any later version. ! ! HFODD is distributed in the hope that it will be useful, but ! WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License ! along with HFODD. If not, see . ! ! OUR NOTICE AND TERMS AND CONDITIONS OF THE GNU GENERAL PUBLIC ! LICENSE ! ! Our Preamble Notice ! ! A. This notice is required to be provided under our contract ! with the U.S. Department of Energy (DOE). This work was ! produced at the Lawrence Livermore National Laboratory under ! Contract No. DE-AC52-07NA27344 with the DOE. ! B. Neither the United States Government nor Lawrence Livermore ! National Security, LLC nor any of their employees, makes any ! warranty, express or implied, or assumes any liability or ! responsibility for the accuracy, completeness, or usefulness ! of any information, apparatus, product, or process disclosed, ! or represents that its use would not infringe privately-owned ! rights. ! C. Also, reference herein to any specific commercial products, ! process, or services by trade name, trademark, manufacturer ! or otherwise does not necessarily constitute or imply its ! endorsement, recommendation, or favoring by the United States ! Government or Lawrence Livermore National Security, LLC. The ! views and opinions of authors expressed herein do not ! necessarily state or reflect those of the United States ! Government or Lawrence Livermore National Security, LLC, and ! shall not be used for advertising or product endorsement ! purposes. ! ! The precise terms and conditions for copying, distribution and ! modification are contained in the file COPYING. ! !*********************************************************************** !----------------------------------------------------------------------! ! ! ! MPI Management Environment for HFODD ! ! ! ! Authors: N. Schunck, ORNL ! ! ! ! This module is/will be in charge of: ! ! - defining the list of tasks to be run based on sequential and ! ! parallel inputs read in hfodd_mpiio.f90 ! ! - managing the load balancing of the tasks ! ! ! ! First included in official release v249e ! ! ! !----------------------------------------------------------------------! Module hfodd_mpimanager Use hfodd_mpiio Implicit None ! Public Integer :: NOFJOB ! Private Logical :: filled Integer :: debug = 0 Integer :: numberQ=4,numberQ_restart=4,numberQ_orig=4, & numberConstraints_orig=1,neck_numberPoints_orig=1 Integer, Allocatable :: IZ_JOB(:),IN_JOB(:),indexQjob(:),indexForce(:),indexTemperature(:), & indexRestartJob(:),indexNeck(:) Integer, Allocatable :: indexSizeJob(:),indexStatesJob(:),indexLengthJob(:),idxDeformJob(:), & basis_size(:),basis_states(:),flex_lambda(:),flex_miu(:),index_file(:), & mapping_qlm(:),mapping_qlm_new(:),mapping_neck(:) Integer, Allocatable :: lambdaJob(:,:),miuJob(:,:),lambdaJob_restart(:,:),miuJob_restart(:,:) Integer, Allocatable :: table_fichiers(:) Real :: neck_step=0.0,neck_step_restart Real, Allocatable :: XIZJOB(:),XINJOB(:) Real, Allocatable :: flex_qlm(:),basis_freq(:),basis_defo(:), & neck_value(:),neck_value_restart(:) Real, Allocatable :: qValueJob(:,:),qValueJob_restart(:,:) Character(Len=68), Allocatable :: fichier(:),fichier_tho(:) Character(Len=68), Allocatable :: fichier_rec_old(:),fichier_rec_new(:),& fichier_tho_old(:),fichier_tho_new(:),& fichier_lic_old(:),fichier_lic_new(:) Character(Len=68), Allocatable :: fichier_rec_tmp(:),& fichier_tho_tmp(:),& fichier_lic_tmp(:) Contains !---------------------------------------------------------------------! ! Defines default deformation grid for 0<= modeMPI_deformation <=2 ! ! Outputs: - numberQ, total number of points on the mesh ! ! - lambdaJob(:), miuJob(:), qValueJob(:), vectors storing ! ! the lambda, miu and q(lambda,miu) values at each point ! ! of the mesh ! ! - numberQ_orig, numberConstraints_orig store the sizes of ! ! the mesh (total size and number of constraints, resp.) ! ! for use when modeMPI_deformation = 2 ! !---------------------------------------------------------------------! Subroutine define_gridQlm() Integer :: iq,nQtotal,nQfactor,ii,lambda,miu,nQ,mpi_rank,mpi_err Integer, Allocatable :: nQmodul(:) Real, Allocatable :: deltaQ(:) Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err) ! Determine total number of points on the mesh: number of collective ! variables 'numberConstraints' times the number of points for each numberQ = 1 Do iq=1,numberConstraints numberQ = numberQ*qStepNumber(iq) End Do If(debug.Ge.1 .And. mpi_rank.Eq.0) & Write(6,'("numberQ = ",i10," numberConstraints=",i10)') numberQ,numberConstraints ! Storing characteristics of the default mesh for future use (if ! modeMPI_deformation = 2) numberQ_orig = numberQ numberConstraints_orig = numberConstraints ! Getting the meshsize for each collective variable Allocate(deltaQ(1:numberConstraints),nQmodul(0:numberConstraints)) Do iq=1,numberConstraints If(qstepNumber(iq) > 1 ) Then deltaQ(iq) = (qFinValue(iq) - qinitValue(iq)) / Real(qstepNumber(iq) - 1) Else deltaQ(iq) = 0.0d0 End If End Do ! Defining the mesh: lambdaJob(ii,iq), miuJob(ii,iq) and qValueJob(ii,iq) give, ! respectively, the value of lambda, miu and q(lambda,miu) at point ii of the ! total mesh and for the collective variable iq Allocate(qValueJob(1:numberQ,1:numberConstraints)) Allocate(lambdaJob(1:numberQ,1:numberConstraints),miuJob(1:numberQ,1:numberConstraints)) nQtotal=1; nQmodul(0)=1; nQfactor=0 Do iq=1,numberConstraints nQ=qStepNumber(iq) lambda=Abs(qLambda(iq)); miu=qMiu(iq) nQtotal=nQtotal*nQ; nQmodul(iq)=numberQ/nQtotal Do ii=1,numberQ lambdaJob(ii,iq)=lambda miuJob(ii,iq)=miu qValueJob(ii,iq)=qinitValue(iq) & + deltaQ(iq) * ((ii-1)/nQmodul(iq) - nQfactor*((ii-1)/nQmodul(iq-1))*nQ) End Do nQfactor = 1 End Do Deallocate(deltaQ,nQmodul) End Subroutine define_gridQlm !---------------------------------------------------------------------! ! Defines restart deformation grid for modeMPI_deformation = 2, 3 ! ! Outputs: - numberQ_restart, total number of points on the restart ! ! mesh ! ! - lambdaJob_restart(:), miuJob_restart(:), and ! ! qValueJob_restart(:), vectors storing the lambda, miu ! ! and q(lambda,miu) values at each point of the restart ! ! mesh ! !---------------------------------------------------------------------! Subroutine define_gridQlm_new() Integer :: iq,nQtotal,nQfactor,ii,lambda,miu,nQ,mpi_rank,mpi_err Integer, Allocatable :: nQmodul(:) Real, Allocatable :: deltaQ(:) Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err) ! Determine total number of points on the restart mesh numberQ_restart = 1 Do iq=1,numberConstraints_restart numberQ_restart = numberQ_restart * qStepNumber_restart(iq) End Do If(debug.Ge.1 .And. mpi_rank.Eq.0) & Write(6,'("numberQ_restart = ",i10," numberConstraints_restart=",i10)') & numberQ_restart,numberConstraints_restart ! Getting the meshsize for each collective variable on the restart mesh Allocate(deltaQ(1:numberConstraints_restart),nQmodul(0:numberConstraints_restart)) Do iq=1,numberConstraints_restart If(qstepNumber_restart(iq)>1) Then deltaQ(iq) = (qFinValue_restart(iq) - qinitValue_restart(iq)) / Real(qstepNumber_restart(iq) - 1) Else deltaQ(iq) = 0.0D0 End If End Do ! Defining the restart mesh, see subroutine define_gridQlm for explanations Allocate(qValueJob_restart(1:numberQ_restart,1:numberConstraints_restart)) Allocate(lambdaJob_restart(1:numberQ_restart,1:numberConstraints_restart),& miuJob_restart(1:numberQ_restart,1:numberConstraints_restart)) nQtotal=1; nQmodul(0)=1; nQfactor=0 Do iq=1,numberConstraints_restart nQ=qStepNumber_restart(iq) lambda=Abs(qLambda_restart(iq)); miu=qMiu_restart(iq) nQtotal=nQtotal*nQ; nQmodul(iq)=numberQ_restart / nQtotal Do ii=1,numberQ_restart lambdaJob_restart(ii,iq)=lambda miuJob_restart(ii,iq)=miu qValueJob_restart(ii,iq)=qinitValue_restart(iq) & + deltaQ(iq) * ( (ii-1)/nQmodul(iq) - nQfactor*((ii-1)/nQmodul(iq-1))*nQ ) End Do nQfactor = 1 End Do Deallocate(deltaQ,nQmodul) End Subroutine define_gridQlm_new !---------------------------------------------------------------------! ! Defines the deformation path for modeMPI_deformation = 4. The full ! ! deformation path is stored in the file hfodd_path.d. It is read in ! ! this routine by processor 0 and broadcast to all processors. ! ! Outputs: - numberQ_restart, total number of points in the path ! ! - lambdaJob_restart(:), miuJob_restart(:), and ! ! qValueJob_restart(:), vectors storing the lambda, miu ! ! and q(lambda,miu) values at each point along the path ! !---------------------------------------------------------------------! Subroutine define_pathQlm_new() Integer, Parameter :: ndatin=84 Integer :: i,lines,iline,number_cons,ierr,mpi_rank,mpi_err Character(Len=80) :: filein ! Getting the rank of the current process Call mpi_comm_rank(MPI_COMM_WORLD, mpi_rank, mpi_err) If(mpi_rank.Eq.0) Then ! open the file filein = "hfodd_path.d" Open(ndatin,file=filein,status='old',form='formatted',iostat=ierr) If(ierr.Ne.0) Then Write(6,'("Error in opening file : ",a80)') filein Stop 'Error in define_pathQlm_new - I/O' End if ! Read the number of constraints used to define a trajectory Read(ndatin,*) number_cons, lines If(debug.Ge.1) Write(6,'("Number of constraints: ",i10," Number of lines: ",i10)') number_cons,lines numberConstraints_restart=number_cons numberQ_restart=lines ! Read trajectory in N-d space (N=number_cons) Allocate(qValueJob_restart(1:numberQ_restart,1:numberConstraints_restart)) Allocate(lambdaJob_restart(1:numberQ_restart,1:numberConstraints_restart),& miuJob_restart(1:numberQ_restart,1:numberConstraints_restart)) iline=0 Do While(iline < lines) iline=iline+1 Read(ndatin,*) (lambdaJob_restart(iline,i),miuJob_restart(iline,i),qValueJob_restart(iline,i),& i=1,number_cons) End Do ! Close file Close(ndatin) End If ! Broadcast number of constraints and number of lines Call MPI_Bcast(numberConstraints_restart,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpi_err) Call MPI_Bcast(numberQ_restart,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpi_err) ! Allocate space for vectors If(mpi_rank.Ne.0) Then Allocate(qValueJob_restart(1:numberQ_restart,1:numberConstraints_restart)) Allocate(lambdaJob_restart(1:numberQ_restart,1:numberConstraints_restart),& miuJob_restart(1:numberQ_restart,1:numberConstraints_restart)) End If ! Broadcast lambda values , Call MPI_Bcast(lambdaJob_restart,numberConstraints_restart*numberQ_restart,MPI_INTEGER,0,& MPI_COMM_WORLD,mpi_err) ! Broadcast miu values , Call MPI_Bcast(miuJob_restart,numberConstraints_restart*numberQ_restart,MPI_INTEGER,0,& MPI_COMM_WORLD,mpi_err) ! Broadcast qlm values Call MPI_Bcast(qValueJob_restart,numberConstraints_restart*numberQ_restart,MPI_DOUBLE_PRECISION,0,& MPI_COMM_WORLD,mpi_err) If(debug.Ge.1) Write(6,'("In define_pathQLM_new - mpi_rank = ",i10, & & " has received numberConstraints_restart = ",i10)') & mpi_rank, numberConstraints_restart End Subroutine define_pathQlm_new !---------------------------------------------------------------------! ! Defines default grid for the constraint on the size of the neck QN ! ! in the case modeMPI_neck >= 1 ! ! Outputs: - neck_value(:), vector storing the values of QN ! ! - neck_numberPoints_orig stores the number of QN values ! ! for use when modeMPI_neck = 2 ! !---------------------------------------------------------------------! Subroutine define_gridNeck() Integer :: i_neck neck_numberPoints_orig = neck_numberPoints Allocate(neck_value(1:neck_numberPoints)) If(neck_numberPoints>1) Then neck_step = (neck_finValue - neck_initValue) / Real(neck_numberPoints - 1) do i_neck=1,neck_numberPoints neck_value(i_neck) = neck_initValue + Real(i_neck - 1)*neck_step End do else neck_value(1) = neck_initValue End if End Subroutine define_gridNeck !---------------------------------------------------------------------! ! Defines restart grid for the constraint on the size of the neck QN ! ! in the case modeMPI_neck = 2 ! ! Outputs: - neck_value_restart(:), vector storing the values of QN ! ! on the restart grid ! !---------------------------------------------------------------------! Subroutine define_gridNeck_new() Integer :: i_neck Allocate(neck_value_restart(1:neck_numberPoints_restart)) If(neck_numberPoints_restart>1) Then neck_step_restart = (neck_finValue_restart - neck_initValue_restart) & / Real(neck_numberPoints_restart - 1) do i_neck=1,neck_numberPoints_restart neck_value_restart(i_neck) = neck_initValue_restart + Real(i_neck - 1)*neck_step_restart End do else neck_value_restart(1) = neck_initValue_restart End if End Subroutine define_gridNeck_new !---------------------------------------------------------------------! ! Defines default deformation grid for modeMPI_deformation >= 3. In ! ! this mode, the initial mesh is extracted on-the-fly from a list of ! ! HFODD record files. As a result, the mesh can be irregular and ! ! have different number of constraints from the restart mesh. The ! ! reading of the HFODD record files is done in parallel to handle ! ! large file sets (and reducing the overhead of reading HFODD binary ! ! files). ! ! ATTENTION: We make the assumption that all files in the list have ! ! the same number (=numberConstraints) and type (qLambda, ! ! qMiu arrays) of constraints. This information is extrac-! ! ted from the MULTICONS keyword. It is necessary to spe- ! ! cify what the constraints are, even though the actual ! ! mesh will be discarded for these two modes. ! ! Outputs: - numberQ, total number of points on the mesh ! ! - lambdaJob(:), miuJob(:), qValueJob(:), vectors storing ! ! the lambda, miu and q(lambda,miu) values at each point ! ! of the mesh ! ! - numberQ_orig, numberConstraints_orig store the sizes of ! ! the mesh (total size and number of constraints, resp.) ! ! - neck_numberPoints_orig is set to 1, as this mode is not ! ! yet compatible with any constraint on the neck. ! !---------------------------------------------------------------------! Subroutine define_restartPointsQlm(NUMITE,NXHERM,NYHERM,NZHERM,IPCONT,ILCONT,IACONT, & IMCONT,IRCONT,IRENMA,IRENIN,LIPKIN,LIPKIP,REFERN, & REFERP,REDELN,REDELP,REFE2N,REFE2P,IFIBLN,INIBLN, & IFIBLP,INIBLP,NMUMAX,NSIMAX,NMUCON,ISHIFT,IF_THO, & IERROR) Integer, INTENT(IN) :: NUMITE,NXHERM,NYHERM,NZHERM,IPCONT,ILCONT,IACONT,IMCONT, & IRCONT,IRENMA,IRENIN,LIPKIN,LIPKIP,IFIBLN,INIBLN,IFIBLP, & INIBLP,NMUMAX,NSIMAX,NMUCON,ISHIFT,IF_THO Integer, INTENT(INOUT) :: IERROR Real, INTENT(INOUT) :: REFERN,REFERP,REDELN,REDELP,REFE2N,REFE2P Logical :: INCREP Integer :: mpi_rank,mpi_err,mpi_world,status(MPI_STATUS_SIZE) Integer :: NFIREP,tag,counter,numberConstraintsMax,exitstat Integer :: number_orig,task,iq,stringLength=68,size_array,mpi_string_68 Integer, Allocatable :: buf_lambdaJob(:),buf_miuJob(:),l_temp(:,:),m_temp(:,:) Character(Len=6) :: string_number Character(Len=68) :: FILREP Character(Len=150) :: commande Real, Allocatable :: buf_qValueJob(:),q_temp(:,:) Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err) Call mpi_comm_size(MPI_COMM_WORLD,mpi_world,mpi_err) ! Define MPI type for a character string of length stringLength=68 and commit it Call mpi_type_contiguous(stringLength,MPI_CHARACTER,mpi_string_68,mpi_err) Call mpi_type_commit(mpi_string_68,mpi_err) ! Process 0 reads list of files to use as restart If(mpi_rank.Eq.0) Then numberQ=1; Call read_fileList() If(debug.Ge.1 .And. mpi_rank.Eq.0) Write(6,'("From hfodd_files.d - numberQ = ",i10)') numberQ End If ! Broadcasting size of collective space, and allocate space for list of filenames (at this ! point, only process 0 has done the allocation) Call MPI_Bcast(numberQ,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpi_err) If(mpi_rank.Ne.0) Then Allocate(fichier(1:numberQ),index_file(1:numberQ)) End If Call mpi_barrier(MPI_COMM_WORLD,mpi_err) ! Broadcasting list of files Call MPI_Bcast(fichier,numberQ,mpi_string_68,0,MPI_COMM_WORLD,mpi_err) Call MPI_Bcast(index_file,numberQ,MPI_INTEGER,0,MPI_COMM_WORLD,mpi_err) ! Allocate space for THO files if needed. Assume there is a directory ./tho/ ! with all required files (this is the consequence of being in restart mode) If(IF_THO.Ge.1) Then Allocate(fichier_tho(1:numberQ)) Do counter=1,numberQ Write(string_number,'(i6.6)') index_file(counter) fichier_tho(counter)='./tho/'//string_number//'.hel' End Do End If ! Make a backup copy of the requested file somewhere else Do counter=1,numberQ If(Mod(mpi_rank,mpi_world).Eq.Mod(counter-1,mpi_world)) Then commande = 'cp '//Trim(fichier(counter))//' .'; exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," counter = ",i10," - Executed command: ",a)') & mpi_rank,counter,commande End If End Do Call mpi_barrier(MPI_COMM_WORLD,mpi_err) ! Reading files and obtaining characteristics of restart jobs. MPI tags are used to synchronize ! the order of send and receives If(mpi_rank.Eq.0) Then Allocate(qValueJob(1:numberQ,1:numberConstraints)) Allocate(lambdaJob(1:numberQ,1:numberConstraints),miuJob(1:numberQ,1:numberConstraints)) qValueJob(:,:)=0.0D0; lambdaJob(:,:)=0; miuJob(:,:)=0 ! Receive the multipole constraints read from each file Do task=1,numberQ tag=task ! Allocate space for current buffers Allocate(buf_lambdaJob(numberConstraints),& buf_miuJob(numberConstraints), & buf_qValueJob(numberConstraints)) If(debug.Ge.3) Write(6,'("mpi_rank ",i10," has allocated buf_[]")') mpi_rank ! Receive local buffers tag= numberQ+task; Call mpi_recv(buf_lambdaJob,numberConstraints,MPI_INTEGER,MPI_ANY_SOURCE,tag, & MPI_COMM_WORLD,status,mpi_err) tag=2*numberQ+task; Call mpi_recv(buf_miuJob,numberConstraints,MPI_INTEGER,MPI_ANY_SOURCE,tag, & MPI_COMM_WORLD,status,mpi_err) tag=3*numberQ+task; Call mpi_recv(buf_qValueJob,numberConstraints,MPI_DOUBLE_PRECISION,MPI_ANY_SOURCE, & tag,MPI_COMM_WORLD,status,mpi_err) If(debug.Ge.3) Write(6,'("mpi_rank ",i10," has received buf_[]")') mpi_rank lambdaJob(task,1:numberConstraints)=buf_lambdaJob(1:numberConstraints) miuJob(task,1:numberConstraints)=buf_miuJob(1:numberConstraints) qValueJob(task,1:numberConstraints)=buf_qValueJob(1:numberConstraints) Deallocate(buf_lambdaJob,buf_miuJob,buf_qValueJob) End Do Else Do task=1,numberQ If(Mod(mpi_rank,mpi_world).Eq.Mod(task-1,mpi_world)) Then tag=task FILREP=fichier(task); NFIREP=42; INCREP=.False. If(debug.Ge.3) Write(6,'("mpi_rank ",i10,", task ",i10,", index_file = ",i10)') & mpi_rank,task,index_file(task) ! Read each file Call read_hfodd(NFIREP,FILREP,NUMITE,NXHERM,NYHERM,NZHERM,IPCONT, & ILCONT,IACONT,IMCONT,IRCONT,IRENMA,IRENIN,LIPKIN, & LIPKIP,REFERN,REFERP,REDELN,REDELP,REFE2N,REFE2P, & IFIBLN,INIBLN,IFIBLP,INIBLP,NMUMAX,NSIMAX,NMUCON, & ISHIFT,IERROR) ! Make sure to close file if not done before Inquire(unit=NFIREP,opened=INCREP); If(INCREP) Close(NFIREP) ! Allocate space for current buffers Allocate(buf_lambdaJob(numberConstraints),& buf_miuJob(numberConstraints), & buf_qValueJob(numberConstraints)) ! For each constraint iq, flex_lambda(iq) contains the lambda value, flex_miu(iq) the miu ! value and flex_qlm(iq) the q(lambda,miu) value buf_lambdaJob(:)=flex_lambda(:) buf_miuJob(:)=flex_miu(:) buf_qValueJob(:)=flex_qlm(:) tag= numberQ+task; Call mpi_send(buf_lambdaJob,numberConstraints,MPI_INTEGER,0,tag, & MPI_COMM_WORLD,mpi_err) tag=2*numberQ+task; Call mpi_send(buf_miuJob,numberConstraints,MPI_INTEGER,0,tag, & MPI_COMM_WORLD,mpi_err) tag=3*numberQ+task; Call mpi_send(buf_qValueJob,numberConstraints,MPI_DOUBLE_PRECISION,0, & tag,MPI_COMM_WORLD,mpi_err) Deallocate(buf_lambdaJob,buf_miuJob,buf_qValueJob) If(debug.Ge.3) Write(6,'("mpi_rank ",i10," has sent buf_[]")') mpi_rank End If End Do End If If(debug.Ge.3) Write(6,'("mpi_rank = ",i10," has gone over the main loop")') mpi_rank If(.Not.Allocated(qValueJob)) Then Allocate(qValueJob(1:numberQ,1:numberConstraints)) Allocate(lambdaJob(1:numberQ,1:numberConstraints),miuJob(1:numberQ,1:numberConstraints)) End If If(debug.Ge.1) Write(6,'("mpi_rank = ",i10," has numberConstraints = ",i10)') & mpi_rank, numberConstraints Call mpi_barrier(MPI_COMM_WORLD,mpi_err) ! Broadcast the full set of multipole constraints to all. First flatten 2D arrays into large vectors size_array=numberQ*numberConstraints Allocate(buf_lambdaJob(1:size_array),buf_miuJob(1:size_array),buf_qValueJob(1:size_array)) If(mpi_rank.Eq.0) Then counter = 0 Do task=1,numberQ Do iq=1,numberConstraints counter=counter+1 buf_lambdaJob(counter)=lambdaJob(task,iq) buf_miuJob(counter)=miuJob(task,iq) buf_qValueJob(counter)=qValueJob(task,iq) End Do End Do Else buf_lambdaJob(:)=0; buf_miuJob(:)=0; buf_qValueJob(:)=0.0D0 End If Call MPI_Bcast(buf_lambdaJob,size_array,MPI_INTEGER,0,MPI_COMM_WORLD,mpi_err) Call MPI_Bcast(buf_miuJob,size_array,MPI_INTEGER,0,MPI_COMM_WORLD,mpi_err) Call MPI_Bcast(buf_qValueJob,size_array,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,mpi_err) ! Then reconstruct 2D arrays from broadcast vectors counter = 0 Do task=1,numberQ Do iq=1,numberConstraints counter=counter+1 lambdaJob(task,iq)=buf_lambdaJob(counter) miuJob(task,iq)=buf_miuJob(counter) qValueJob(task,iq)=buf_qValueJob(counter) End Do End Do Deallocate(buf_lambdaJob,buf_miuJob,buf_qValueJob) If(debug.Ge.4) Then If(mpi_rank.Eq.0) Then Do task=1,numberQ Write(6,'("File number",i10," is ",a68)') task,fichier(task) Do iq=1,numberConstraints Write(6,'("l = ",i4," m = ",i4," qlm = ",f10.5)') & lambdaJob(task,iq),miuJob(task,iq),qValueJob(task,iq) End Do End Do End If End If ! Set original number of constraints and original number of "grid" points numberQ_orig = numberQ numberConstraints_orig = numberConstraints neck_numberPoints = 1 ! Releasing derived MPI datatype Call mpi_type_free(mpi_string_68,mpi_err) End Subroutine define_restartPointsQlm !---------------------------------------------------------------------! ! Function returning the job number embedded in a HFODD record file ! ! name, used when modeMPI_deformation = 3 only ! !---------------------------------------------------------------------! Integer Function extract_index(filename) Character(Len=68), INTENT(IN) :: filename Integer :: i If(filename(1:1).Eq.'H') Then Read(filename,'("HFODD_",i8.8,".REC")') i extract_index = i End If End Function extract_index !---------------------------------------------------------------------! ! Subroutine reading the list of files to restart from in the case ! ! where modeMPI_deformation = 3. Executed only from proc 0 ! ! Outputs: - numberQ, total number of points on the mesh ! ! - fichier(:), list of files to restart from ! !---------------------------------------------------------------------! Subroutine read_fileList() Integer :: unit,counter Character(Len=68) :: filename unit=41 ! First pass through the file to count the number of files: defines the initial "mesh" counter=0; filename='' Open(unit,FILE='hfodd_files.d',STATUS='UNKNOWN',FORM='FORMATTED') Do While(trim(filename) .Ne. 'ALL_DONE' .And. trim(filename) .Ne. 'ALL-DONE') Read(unit,*) filename counter=counter+1 End Do Close(unit) numberQ=counter-1 ! Second pass to actually read the filename filename='' Allocate(fichier(1:numberQ),index_file(1:numberQ)) Open(unit,FILE='hfodd_files.d',STATUS='UNKNOWN',FORM='FORMATTED') Do counter=1,numberQ Read(unit,*) filename index_file(counter)=extract_index(filename) fichier(counter)='./restart/'//filename(1:60) End Do Close(unit) End Subroutine read_fileList !---------------------------------------------------------------------! ! Subroutine mapping the default deformation grid to the restart ! ! deformation grid. Only used when modeMPI_deformation = 2. For this ! ! mapping to function, it is assumed that the restart collective ! ! space is equal or larger than the original one: there may be more ! ! constraints, but not less; the number of points for the shared ! ! collective variable may be different. ! ! Outputs: - ii=mapping_qlm(jj) gives the index ii of the original ! ! grid used to restart the calculation on point jj of the ! ! restart grid ! ! - numberQ, numberConstraints are overwritten and contain ! ! the characteristics of the restart collective space ! ! - qValueJob(:), lambdaJob(:) and miuJob(:) are similarly ! ! associated with the restart grid ! ! TODO: THE CURRENT METHOD MAY BE TAXING THE OS WHEN THE NEW GRID IS ! ! LARGER THAN THE OLD ONE. INDEED, THE SAME FILE MAY BE READ ! ! BY MANY DIFFERENT PROCESSES, WHICH IS NOT EFFICIENT. ! !---------------------------------------------------------------------! Subroutine mapping_gridQlm(IF_THO,LIPKIN,LIPKIP) Integer, INTENT(IN) :: IF_THO,LIPKIN,LIPKIP Integer :: ii,jj,index_optimal,iq,jq,shift,lambda_restart,miu_restart, & lambda,miu,mpi_rank,mpi_err,mpi_world,exitstat,pass_number Real :: distance_min,distance,qlm_restart,qlm Character(Len=132) :: directory_name Character(Len=68) :: fichier_courant,fichier_optimal Character(Len=150) :: commande Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err) Call mpi_comm_size(MPI_COMM_WORLD,mpi_world,mpi_err) ! Default mapping: all new calculations restart from HFODD_00000001.REC Allocate(mapping_qlm(1:numberQ_restart)) Do jj=1,numberQ_restart mapping_qlm(jj)=1 End Do ! Move all restart files into current directory If(mpi_rank.Eq.0.And.modeMPI_neck.Ne.2) Then commande = 'mv restart restart_old'; exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') & mpi_rank,commande commande = 'mkdir restart'; exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') & mpi_rank,commande If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then commande = 'mv lic lic_old'; exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') & mpi_rank,commande commande = 'mkdir lic'; exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') & mpi_rank,commande End If ! In batch mode, restart files will later be overwritten with record files. We make an ! additional copy here If(batch_mode.Eq.1) Then commande = 'mkdir restart_tmp'; exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') & mpi_rank,commande If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then commande = 'mkdir lic_tmp'; exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') & mpi_rank,commande End If End If End If ! Wait here for all processes Call mpi_barrier(MPI_COMM_WORLD,mpi_err) Allocate(fichier_rec_old(1:numberQ_restart),fichier_rec_new(1:numberQ_restart)) If(IF_THO.Ge.1) Allocate(fichier_tho_old(1:numberQ_restart),fichier_tho_new(1:numberQ_restart)) If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Allocate(fichier_lic_old(1:numberQ_restart),fichier_lic_new(1:numberQ_restart)) If(batch_mode.Eq.1) Then Allocate(fichier_rec_tmp(1:numberQ_restart)) If(IF_THO.Ge.1) Allocate(fichier_tho_tmp(1:numberQ_restart)) If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Allocate(fichier_lic_tmp(1:numberQ_restart)) End If ! Loop over the new grid Do jj=1,numberQ_restart index_optimal=1; distance_min=1.D9 ! Loop over the old grid Do ii=1,numberQ distance = 0.0D0 ! Loop over the old constraints Do iq=1,numberConstraints ! Both old and new constraints were read using linked lists. If Nold is the ! number of old constraints and Nnew that of the new constraints, we have ! Nnew >= Nold (necessarily) and the shared Nold constraints are stored at ! the end of lambdaJob(:), lambdaJob_restart(:), etc. To ge the proper ! correspondence, we use the shift below shift = iq + (numberConstraints_restart - numberConstraints) lambda_restart=lambdaJob_restart(jj,shift) miu_restart= miuJob_restart(jj,shift) qlm_restart=qValueJob_restart(jj,shift) lambda=lambdaJob(ii,iq) miu=miuJob(ii,iq) qlm=qValueJob(ii,iq) ! Measuring the distance between point ii (old) and jj (new) ! in terms of multipole moments If(lambda.Eq.lambda_restart .And. miu.Eq.miu_restart) Then distance = distance + (qlm - qlm_restart)**2 End If End Do ! Looking for the point index_optimal that minimizes the distance If(distance.Le.distance_min) Then distance_min=distance; index_optimal = ii End If End Do ! mv proper file into restart directory with proper name If(modeMPI_neck.Ne.2)Then Write(fichier_optimal,'("./restart_old/HFODD_",I8.8,".REC")') index_optimal Write(fichier_courant,'("./restart/HFODD_",I8.8,".REC")') jj fichier_rec_old(jj)=Trim(fichier_optimal) fichier_rec_new(jj)=Trim(fichier_courant) If(IF_THO.Ge.1) Then Write(fichier_optimal,'("./tho/t",I6.6,".hel")') index_optimal Write(fichier_courant,'("./t",I6.6,".hel")') jj fichier_tho_old(jj)=Trim(fichier_optimal) fichier_tho_new(jj)=Trim(fichier_courant) End If If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then Write(fichier_optimal,'("./lic_old/HFODD_",I8.8,".LIC")') index_optimal Write(fichier_courant,'("./lic/HFODD_",I8.8,".LIC")') jj fichier_lic_old(jj)=Trim(fichier_optimal) fichier_lic_new(jj)=Trim(fichier_courant) End If mapping_qlm(jj)=jj ! In batch mode, we need to make an extra copy of all these restart files, since they will ! be overwritten later If(batch_mode.Eq.1) Then Write(fichier_courant,'("./restart_tmp/HFODD_",I8.8,".REC")') jj fichier_rec_tmp(jj)=Trim(fichier_courant) If(IF_THO.Eq.1) Then Write(fichier_courant,'("./t",I6.6,".hel")') jj fichier_tho_tmp(jj)=Trim(fichier_courant) End If If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then Write(fichier_courant,'("./lic_tmp/HFODD_",I8.8,".LIC")') jj fichier_lic_tmp(jj)=Trim(fichier_courant) End If End If Else mapping_qlm(jj)=index_optimal End If If(mpi_rank.Eq.0.And.debug.Ge.4) Then Write(6,'("mpi_rank = ",i10," mapping_qlm - jj = ",i10," index_optimal = ",i10)') & mpi_rank,jj,index_optimal Do jq=1,numberConstraints_restart Write(6,'("RESTART - mpi_rank = ",i10," l=",i2,"m=",i2," Qlm=",f10.5)') & mpi_rank,lambdaJob_restart(jj,jq), & miuJob_restart(jj,jq), & qValueJob_restart(jj,jq) End do Do iq=1,numberConstraints Write(6,'("ORIGINAL - mpi_rank = ",i10," l=",i2,"m=",i2," Qlm=",f10.5)') & mpi_rank,lambdaJob(index_optimal,iq), & miuJob(index_optimal,iq), & qValueJob(index_optimal,iq) End do End If End Do If(modeMPI_neck.Ne.2) Then ! Loop over the new grid Do jj=1,numberQ_restart If(Mod(mpi_rank,mpi_world).Eq.Mod(jj-1,mpi_world)) Then commande = 'cp '//fichier_rec_old(jj)//' '//fichier_rec_new(jj); exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') & mpi_rank,jj,commande If(IF_THO.Ge.1) Then commande = 'cp '//fichier_tho_old(jj)//' '//fichier_tho_new(jj); exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') & mpi_rank,jj,commande End If If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then commande = 'cp '//fichier_lic_old(jj)//' '//fichier_lic_new(jj); exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') & mpi_rank,jj,commande End If ! Make extra copy in batch mode If(batch_mode.Eq.1) Then commande = 'cp '//fichier_rec_old(jj)//' '//fichier_rec_tmp(jj); exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') & mpi_rank,jj,commande If(IF_THO.Ge.1) Then commande = 'cp '//fichier_tho_old(jj)//' '//fichier_tho_tmp(jj); exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') & mpi_rank,jj,commande End If If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then commande = 'cp '//fichier_lic_old(jj)//' '//fichier_lic_tmp(jj); exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') & mpi_rank,jj,commande End If End If End If End Do ! Wait here for all processes Call mpi_barrier(MPI_COMM_WORLD,mpi_err) ! Check that all files have truly been copied Allocate(table_fichiers(1:numberQ_restart)); table_fichiers(:)=0 directory_name='restart' pass_number=1; filled=.False. Do While (.Not.filled) If(mpi_rank.Eq.0.And.debug.Ge.2) Write(6,'("Pass=",i5)') pass_number Call check_directory(directory_name,numberQ_restart) Call fill_holes(numberQ_restart,IF_THO,LIPKIN,LIPKIP) pass_number=pass_number+1 If(pass_number.Gt.4) Exit ! avoid infinite loops End Do Deallocate(table_fichiers) ! Check that all files in restart_tmp/ have truly been copied (batch mode) If(batch_mode.Eq.1) Then Allocate(table_fichiers(1:numberQ_restart)); table_fichiers(:)=0 directory_name='restart_tmp' pass_number=1; filled=.False. Do While (.Not.filled) If(mpi_rank.Eq.0.And.debug.Ge.2) Write(6,'("Pass =",i5)') pass_number Call check_directory(directory_name,numberQ_restart) Call fill_holes(numberQ_restart,IF_THO,LIPKIN,LIPKIP) pass_number=pass_number+1 If(pass_number.Gt.4) Exit ! avoid infinite loops End Do Deallocate(table_fichiers) End If End If ! End modeMPI_neck.Ne.2 Deallocate(fichier_rec_old,fichier_rec_new) If(IF_THO.Ge.1) Deallocate(fichier_tho_old,fichier_tho_new) If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Deallocate(fichier_lic_old,fichier_lic_new) If(batch_mode.Eq.1) Then Deallocate(fichier_rec_tmp) If(IF_THO.Ge.1) Deallocate(fichier_tho_tmp) If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Deallocate(fichier_lic_tmp) End If Deallocate(qValueJob,lambdaJob,miuJob) numberQ=numberQ_restart numberConstraints=numberConstraints_restart Allocate(qValueJob(1:numberQ,1:numberConstraints)) Allocate(lambdaJob(1:numberQ,1:numberConstraints),miuJob(1:numberQ,1:numberConstraints)) qValueJob(:,:)=qValueJob_restart(:,:) lambdaJob(:,:)=lambdaJob_restart(:,:); miuJob(:,:)=miuJob_restart(:,:) Deallocate(qValueJob_restart,lambdaJob_restart,miuJob_restart) End Subroutine mapping_gridQlm !---------------------------------------------------------------------! ! Subroutine mapping the default deformation grid to the restart ! ! deformation grid. Only used when modeMPI_deformation >= 3. In this ! ! case, the original collective space is arbitrary compared to the ! ! restart one: it may have more or fewer collective variables, the ! ! mesh can be irregular, etc. This is an option meant to provide ! ! flexibility. ! ! Outputs: - ii=mapping_qlm_new(jj) gives as before the index ii of ! ! the original grid used to restart the calculation on ! ! point jj of the restart grid ! ! - numberQ, numberConstraints are overwritten and contain ! ! the characteristics of the restart collective space ! ! - qValueJob(:), lambdaJob(:) and miuJob(:) are similarly ! ! associated with the restart grid ! !---------------------------------------------------------------------! Subroutine mapping_gridQlm_new(IF_THO,LIPKIN,LIPKIP) Integer, INTENT(IN) :: IF_THO,LIPKIN,LIPKIP Integer :: ii,jj,index_optimal,iq,jq,lambda_restart,miu_restart, & lambda,miu,active,mpi_rank,mpi_err,mpi_world,exitstat,& pass_number Real :: distance_min,distance,qlm_restart,qlm Character(Len=132) :: directory_name Character(Len=68) :: fichier_courant,fichier_optimal Character(Len=150) :: commande Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err) Call mpi_comm_size(MPI_COMM_WORLD,mpi_world,mpi_err) ! Default mapping: all new calculations restart from HFODD_00000001.REC Allocate(mapping_qlm_new(1:numberQ_restart)) Do jj=1,numberQ_restart mapping_qlm_new(jj)=1 End Do ! Move all restart files into current directory If(mpi_rank.Eq.0) Then commande = 'mv restart restart_old'; exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') & mpi_rank,commande commande = 'mkdir restart'; exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') & mpi_rank,commande If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then commande = 'mv lic lic_old'; exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') & mpi_rank,commande commande = 'mkdir lic'; exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') & mpi_rank,commande End If ! In batch mode, restart files will later be overwritten with record files. We make an ! additional copy here If(batch_mode.Eq.1) Then commande = 'mkdir restart_tmp'; exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') & mpi_rank,commande If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then commande = 'mkdir lic_tmp'; exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," executed command: ",a)') & mpi_rank,commande End If End If End If ! Wait here for all processes Call mpi_barrier(MPI_COMM_WORLD,mpi_err) ! Allocate filenames Allocate(fichier_rec_old(1:numberQ_restart),fichier_rec_new(1:numberQ_restart)) If(IF_THO.Ge.1) Allocate(fichier_tho_old(1:numberQ_restart),fichier_tho_new(1:numberQ_restart)) If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Allocate(fichier_lic_old(1:numberQ_restart),fichier_lic_new(1:numberQ_restart)) If(batch_mode.Eq.1) Then Allocate(fichier_rec_tmp(1:numberQ_restart)) If(IF_THO.Ge.1) Allocate(fichier_tho_tmp(1:numberQ_restart)) If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Allocate(fichier_lic_tmp(1:numberQ_restart)) End If ! Loop over the new grid Do jj=1,numberQ_restart index_optimal=1; distance_min=1.D9 ! Loop over the old grid Do ii=1,numberQ distance=0.0D0; active=0 ! Loop over the new collective variables (=constraints) Do jq=1,numberConstraints_restart lambda_restart=lambdaJob_restart(jj,jq) miu_restart=miuJob_restart(jj,jq) qlm_restart=qValueJob_restart(jj,jq) ! Loop over the old collective variables (=constraints) Do iq=1,numberConstraints lambda=lambdaJob(ii,iq) miu=miuJob(ii,iq) qlm=qValueJob(ii,iq) ! Measuring the distance between point ii (old) and jj (new) ! in terms of multipole moments If(lambda.Eq.lambda_restart .And. miu.Eq.miu_restart) Then active=active+1 distance = distance + (qlm - qlm_restart)**2 End If End Do End Do ! Looking for the point index_optimal that minimizes the distance If(distance.Le.distance_min.And.active.Gt.0) Then distance_min=distance; index_optimal = ii End If End Do ! move proper file into restart directory with proper name Write(fichier_optimal,'("./restart_old/HFODD_",I8.8,".REC")') index_file(index_optimal) Write(fichier_courant,'("./restart/HFODD_",I8.8,".REC")') jj fichier_rec_old(jj)=Trim(fichier_optimal) fichier_rec_new(jj)=Trim(fichier_courant) If(IF_THO.Ge.1) Then Write(fichier_optimal,'("./tho/t",I6.6,".hel")') index_file(index_optimal) Write(fichier_courant,'("./t",I6.6,".hel")') jj fichier_tho_old(jj)=Trim(fichier_optimal) fichier_tho_new(jj)=Trim(fichier_courant) End If If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then Write(fichier_optimal,'("./lic_old/HFODD_",I8.8,".LIC")') index_optimal Write(fichier_courant,'("./lic/HFODD_",I8.8,".LIC")') jj fichier_lic_old(jj)=Trim(fichier_optimal) fichier_lic_new(jj)=Trim(fichier_courant) End If ! In batch mode, we need to make an extra copy of all these restart files, since they will ! be overwritten later If(batch_mode.Eq.1) Then Write(fichier_courant,'("./restart_tmp/HFODD_",I8.8,".REC")') jj fichier_rec_tmp(jj)=Trim(fichier_courant) If(IF_THO.Eq.1) Then Write(fichier_courant,'("./t",I6.6,".hel")') jj fichier_tho_tmp(jj)=Trim(fichier_courant) End If If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then Write(fichier_courant,'("./lic_tmp/HFODD_",I8.8,".LIC")') jj fichier_lic_tmp(jj)=Trim(fichier_courant) End If End If mapping_qlm_new(jj)=jj If(mpi_rank.Eq.0.And.debug.Ge.4) Then Write(6,'("mpi_rank = ",i10," mapping_qlm_new - jj = ",i10," index_optimal = ",i10)') & mpi_rank,jj,index_file(index_optimal) Do jq=1,numberConstraints_restart Write(6,'("RESTART - mpi_rank = ",i10," l=",i2,"m=",i2," Qlm=",f10.5)') & mpi_rank,lambdaJob_restart(jj,jq), & miuJob_restart(jj,jq), & qValueJob_restart(jj,jq) End do Do iq=1,numberConstraints Write(6,'("ORIGINAL - mpi_rank = ",i10," l=",i2,"m=",i2," Qlm=",f10.5)') & mpi_rank,lambdaJob(index_optimal,iq), & miuJob(index_optimal,iq), & qValueJob(index_optimal,iq) End do End If End Do ! Loop over the new grid Do jj=1,numberQ_restart If(Mod(mpi_rank,mpi_world).Eq.Mod(jj-1,mpi_world)) Then commande = 'cp '//fichier_rec_old(jj)//' '//fichier_rec_new(jj); exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') & mpi_rank,jj,commande If(IF_THO.Ge.1) Then commande = 'cp '//fichier_tho_old(jj)//' '//fichier_tho_new(jj); exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') & mpi_rank,jj,commande End If If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then commande = 'cp '//fichier_lic_old(jj)//' '//fichier_lic_new(jj); exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') & mpi_rank,jj,commande End If ! Make extra copy in batch mode If(batch_mode.Eq.1) Then commande = 'cp '//fichier_rec_old(jj)//' '//fichier_rec_tmp(jj); exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') & mpi_rank,jj,commande If(IF_THO.Ge.1) Then commande = 'cp '//fichier_tho_old(jj)//' '//fichier_tho_tmp(jj); exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') & mpi_rank,jj,commande End If If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then commande = 'cp '//fichier_lic_old(jj)//' '//fichier_lic_tmp(jj); exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') & mpi_rank,jj,commande End If End If End If End Do ! Wait here for all processes Call mpi_barrier(MPI_COMM_WORLD,mpi_err) ! Check that all files in restart/ have truly been copied Allocate(table_fichiers(1:numberQ_restart)); table_fichiers(:)=0 directory_name='restart' pass_number=1; filled=.False. Do While (.Not.filled) If(mpi_rank.Eq.0.And.debug.Ge.2) Write(6,'("Pass =",i5)') pass_number Call check_directory(directory_name,numberQ_restart) Call fill_holes(numberQ_restart,IF_THO,LIPKIN,LIPKIP) pass_number=pass_number+1 If(pass_number.Gt.4) Exit ! avoid infinite loops End Do Deallocate(table_fichiers) ! Check that all files in restart_tmp/ have truly been copied (batch mode) If(batch_mode.Eq.1) Then Allocate(table_fichiers(1:numberQ_restart)); table_fichiers(:)=0 directory_name='restart_tmp' pass_number=1; filled=.False. Do While (.Not.filled) If(mpi_rank.Eq.0.And.debug.Ge.2) Write(6,'("Pass =",i5)') pass_number Call check_directory(directory_name,numberQ_restart) Call fill_holes(numberQ_restart,IF_THO,LIPKIN,LIPKIP) pass_number=pass_number+1 If(pass_number.Gt.4) Exit ! avoid infinite loops End Do Deallocate(table_fichiers) End If Deallocate(fichier_rec_old,fichier_rec_new) If(IF_THO.Ge.1) Deallocate(fichier_tho_old,fichier_tho_new) If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Deallocate(fichier_lic_old,fichier_lic_new) If(batch_mode.Eq.1) Then Deallocate(fichier_rec_tmp) If(IF_THO.Ge.1) Deallocate(fichier_tho_tmp) If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Deallocate(fichier_lic_tmp) End If Deallocate(qValueJob,lambdaJob,miuJob) numberQ=numberQ_restart numberConstraints=numberConstraints_restart If(mpi_rank.Eq.0.And.debug.Ge.2) & Write(6,'("In mapping_gridQlm_new - numberQ =",i10," numberConstraints = ",i5)') numberQ, numberConstraints Allocate(qValueJob(1:numberQ,1:numberConstraints)) Allocate(lambdaJob(1:numberQ,1:numberConstraints),miuJob(1:numberQ,1:numberConstraints)) qValueJob(:,:)=qValueJob_restart(:,:) lambdaJob(:,:)=lambdaJob_restart(:,:); miuJob(:,:)=miuJob_restart(:,:) Deallocate(qValueJob_restart,lambdaJob_restart,miuJob_restart) End Subroutine mapping_gridQlm_new !---------------------------------------------------------------------! ! Subroutine mapping the default neck grid to the restart neck grid. ! ! Outputs: - ii=mapping_neck(jj) gives the index ii of the file used ! ! to restart the calculation of point jj on the restart ! ! grid. ! ! - neck_value(:), neck_numberPoints are overwritten and ! ! now correspond to the restart grid ! !---------------------------------------------------------------------! Subroutine mapping_gridNeck() Integer :: ii,jj,index_optimal Real :: qn,qn_restart,distance,distance_min ! Default mapping: all new calculations restart from HFODD_00000001.REC Allocate(mapping_neck(1:neck_numberPoints_restart)) Do jj=1,neck_numberPoints_restart mapping_neck(jj)=1 End Do Do jj=1,neck_numberPoints_restart index_optimal=1; distance_min=1.0D9 Do ii=1,neck_numberPoints qn_restart = neck_value_restart(jj) qn = neck_value(ii) distance = (qn - qn_restart)**2 If(distance.Le.distance_min) Then distance_min=distance; index_optimal = ii End If End Do mapping_neck(jj)=index_optimal If(debug.Ge.1) Write(6,'("mapping_gridNeck - jj = ",i10," index_optimal = ",i10)') & jj,index_optimal End Do neck_numberPoints = neck_numberPoints_restart Deallocate(neck_value); Allocate(neck_value(1:neck_numberPoints)) neck_value(:) = neck_value_restart(:) Deallocate(neck_value_restart) End Subroutine mapping_gridNeck !---------------------------------------------------------------------! ! Subroutine defining the tasks ! ! Possible combinations: ! ! modeMPI_deformation 1 1 2 2 2 3 ! ! modeMPI_neck 0 1 0 1 2 0 ! !---------------------------------------------------------------------! Subroutine mpi_constructJobList(NUMITE,NXHERM,NYHERM,NZHERM, & IPCONT,ILCONT,IACONT,IMCONT,IRCONT, & IRENMA,IRENIN, & LIPKIN,LIPKIP, & REFERN,REFERP,REDELN, & REDELP,REFE2N,REFE2P, & IFIBLN,INIBLN,IFIBLP,INIBLP, & NMUMAX,NSIMAX,NMUCON,ISHIFT, & IF_THO, & IERROR) Integer, INTENT(IN) :: NUMITE,NXHERM,NYHERM,NZHERM,IPCONT,ILCONT,IACONT, & IMCONT,IRCONT,IRENMA,IRENIN,LIPKIN,LIPKIP,IFIBLN, & INIBLN,IFIBLP,INIBLP,NMUMAX,NSIMAX,NMUCON,ISHIFT, & IF_THO Integer, INTENT(INOUT) :: IERROR Real, INTENT(INOUT) :: REFERN,REFERP,REDELN,REDELP,REFE2N,REFE2P Integer :: mpi_rank, mpi_err Integer :: iz,in,iq,nCommon,nGridOrig,indexJob_orig,iCommon,force,temperature,idxNeck, & idxDefo,iGrid,indexJob,zNumber,nNumber,i_neck,ii,i_shel,i_stat,i_freq,i_defo Integer, Allocatable :: zVector(:),nVector(:),mapping_oldGrid(:,:),indexJobVector_orig(:,:) Real :: xzNumber, xnNumber Real, Allocatable :: xzVector(:),xnVector(:) Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err) If(debug.Ge.1 .And. mpi_rank .Eq.0) & Write(6,'("modeMPI_deformation=",i3," modeMPI_neck=",i3," modeMPI_basis=",i3)') & modeMPI_deformation,modeMPI_neck,modeMPI_basis ! Define deformation and neck grids If(modeMPI_deformation.Le.2) Then Call define_gridQlm() If(modeMPI_deformation.Eq.2) Then Call define_gridQlm_new() Call mapping_gridQlm(IF_THO,LIPKIN,LIPKIP) End If ! If constraints on the size of the neck are active, we define the grids If(modeMPI_neck.Ge.1) Then Call define_gridNeck() If(modeMPI_neck.Eq.1) Then ! Default mapping if no constraint on the size of the neck is requested, neck_numberPoints ! is preset to 1 in routine mpi_initializeParallelData() of module hfodd_mpiio.f90 neck_numberPoints_orig=1 Allocate(mapping_neck(1:neck_numberPoints)) Do ii=1,neck_numberPoints mapping_neck(ii)=1 End Do Else Call define_gridNeck_new() Call mapping_gridNeck() End If Else ! Default mapping if no constraint on the size of the neck is requested, neck_numberPoints ! is preset to 1 in routine mpi_initializeParallelData() of module hfodd_mpiio.f90 neck_numberPoints_orig=1 Allocate(mapping_neck(1:neck_numberPoints)) Do ii=1,neck_numberPoints mapping_neck(ii)=1 End Do End If Else Call define_restartPointsQlm(NUMITE,NXHERM,NYHERM,NZHERM,IPCONT,ILCONT,IACONT, & IMCONT,IRCONT,IRENMA,IRENIN,LIPKIN,LIPKIP,REFERN, & REFERP,REDELN,REDELP,REFE2N,REFE2P,IFIBLN,INIBLN, & IFIBLP,INIBLP,NMUMAX,NSIMAX,NMUCON,ISHIFT,IF_THO, & IERROR) If(modeMPI_deformation.Eq.3) Call define_gridQlm_new() If(modeMPI_deformation.Eq.4) Call define_pathQlm_new() Call mapping_gridQlm_new(IF_THO,LIPKIN,LIPKIP) ! Default mapping since no constraint on the size of the neck is allowed neck_numberPoints_orig=neck_numberPoints Allocate(mapping_neck(1:neck_numberPoints)) Do ii=1,neck_numberPoints mapping_neck(ii)=ii End Do End If ! default allocations of mapping arrays: if no restart is required, these mappings simply ! provide identity vectors (v(i) = i) If(modeMPI_deformation.Le.1) Then Allocate(mapping_qlm(1:numberQ)) Do ii=1,numberQ mapping_qlm(ii)=ii If(debug.Ge.4 .And. mpi_rank .Eq.0) Write(6,'("ii = ",i10," mapping = ",i10)') ii,mapping_qlm(ii) End Do End If ! Define nucleus grid: integers Allocate(zVector(1:proton_stepNumber),nVector(1:neutron_stepNumber)) Do iz=1,proton_stepNumber zVector(iz) = proton_initValue + (iz-1)*proton_stepValue End Do Do in=1, neutron_stepNumber nVector(in) = neutron_initValue + (in-1)*neutron_stepValue End Do ! Define nucleus grid: float If(x_proton_initValue.Gt.0.0D0) Then Allocate(xzVector(1:proton_stepNumber)) Do iz=1,proton_stepNumber xzVector(iz) = x_proton_initValue + (iz-1)*x_proton_stepValue End Do End If If(x_neutron_initValue.Gt.0.0D0) Then Allocate(xnVector(1:neutron_stepNumber)) Do in=1, neutron_stepNumber xnVector(in) = x_neutron_initValue + (in-1)*x_neutron_stepValue End Do End If ! Default deformation mode If(modeMPI_deformation.Ge.1.And.modeMPI_deformation.Le.2) Then NOFJOB = proton_stepNumber * neutron_stepNumber * counterForce * numberQ & * counterTemperature * neck_numberPoints Allocate(IZ_JOB(1:NOFJOB),IN_JOB(1:NOFJOB),indexQjob(1:NOFJOB),indexRestartJob(1:NOFJOB),& indexForce(1:NOFJOB),indexTemperature(1:NOFJOB),indexNeck(1:NOFJOB)) Allocate(XIZJOB(1:NOFJOB),XINJOB(1:NOFJOB)); XIZJOB(:)=-1.0D0; XINJOB(:)=-1.0D0 nCommon = proton_stepNumber * neutron_stepNumber * counterForce * counterTemperature nGridOrig = numberQ * neck_numberPoints_orig Allocate(indexJobVector_orig(1:nGridOrig, 1:nCommon)) indexJob_orig = 0 If(debug.Ge.1 .And. mpi_rank .Eq.0) Then Write(6,'("Sanity checks:")') Write(6,'("proton_stepNumber=",i10," neutron_stepNumber=",i10," numberConstraints=",i10)') & proton_stepNumber,neutron_stepNumber,numberConstraints Write(6,'("counterForce=",i10," counterTemperature=",i10," restartTemperature=",i10)') & counterForce,counterTemperature,restartTemperature Write(6,'("numberQ=",i10," numberQ_orig=",i10," neck_numberPoints=",i10," neck_numberPoints_orig=",i10)') & numberQ,numberQ_orig,neck_numberPoints,neck_numberPoints_orig End If If(restartTemperature == 0) Then iCommon=0 Allocate(mapping_oldGrid(1:numberQ,1:neck_numberPoints_orig)) Do iz=1,proton_stepNumber Do in=1,neutron_stepNumber Do force=1,counterForce Do temperature=1,counterTemperature iCommon=iCommon+1; iGrid=0 Do idxNeck=1,neck_numberPoints_orig Do idxDefo=1,numberQ iGrid=iGrid+1 mapping_oldGrid(idxDefo,idxNeck)=iGrid indexJob_orig=indexJob_orig+1 indexJobVector_orig(iGrid,iCommon)=indexJob_orig End Do End Do End Do End Do End Do End Do iCommon=0; indexJob=0 Do iz=1,proton_stepNumber zNumber = zVector(iz) If(x_proton_initValue.Gt.0.0D0) xzNumber=xzVector(iz) Do in=1,neutron_stepNumber nNumber = nVector(in) If(x_neutron_initValue.Gt.0.0D0) xnNumber=xnVector(in) Do force=1,counterForce iCommon = iCommon + 1 Do temperature=1,counterTemperature Do i_neck=1,neck_numberPoints Do ii=1,numberQ indexJob=indexJob+1 indexQjob(indexJob)=ii indexForce(indexJob)=force indexTemperature(indexJob)=temperature indexNeck(indexJob)=i_neck IZ_JOB(indexJob)=zNumber IN_JOB(indexJob)=nNumber If(x_proton_initValue.Gt.0.0D0) Then XIZJOB(indexJob)=xzNumber Else XIZJOB(indexJob)=Real(zNumber,Kind=Kind(1.0)) End If If(x_neutron_initValue.Gt.0.0D0) Then XINJOB(indexJob)=xnNumber Else XINJOB(indexJob)=Real(nNumber,Kind=Kind(1.0)) End If idxDefo=mapping_qlm(ii) idxNeck=mapping_neck(i_neck) iGrid=mapping_oldGrid(idxDefo,idxNeck) indexRestartJob(indexJob)=indexJobVector_orig(iGrid,iCommon) If(debug.Ge.4 .And. mpi_rank .Eq.0) Then Do iq=1,numberConstraints Write(6,'("Z=",i3," N=",i3," force=",a4," T=",f8.4, & & " iq=",i2," l=",i2, & & " m=",i2," Qlm=",f10.5," -- RESTART i=",i6)') & zNumber,nNumber, & forceVector(indexForce(indexJob)), & temperatureVector(temperature), & iq,lambdaJob(indexQjob(indexJob),iq), & miuJob(indexQjob(indexJob),iq), & qValueJob(indexQjob(indexJob),iq), & indexJobVector_orig(iGrid,iCommon) End Do End If End Do End Do End Do End Do End Do End Do Deallocate(indexJobVector_orig,mapping_oldGrid,mapping_qlm,mapping_neck) Else iCommon = 0 Allocate(mapping_oldGrid(1:numberQ,1:neck_numberPoints_orig)) Do iz=1,proton_stepNumber Do in=1,neutron_stepNumber Do force=1,counterForce Do temperature=1,counterTemperature iCommon=iCommon+1; iGrid=0 Do idxNeck=1,neck_numberPoints_orig Do idxDefo=1,numberQ iGrid=iGrid+1 mapping_oldGrid(idxDefo,idxNeck)=iGrid indexJob_orig=indexJob_orig+1 indexJobVector_orig(iGrid,iCommon)=indexJob_orig End Do End Do End Do End Do End Do End Do iCommon=0; indexJob=0 Do iz=1,proton_stepNumber zNumber = zVector(iz) If(x_proton_initValue.Gt.0.0D0) xzNumber=xzVector(iz) Do in=1,neutron_stepNumber nNumber = nVector(in) If(x_neutron_initValue.Gt.0.0D0) xnNumber=xnVector(in) Do force=1,counterForce Do temperature=1,counterTemperature iCommon = iCommon + 1 Do i_neck=1,neck_numberPoints Do ii=1,numberQ indexJob=indexJob+1 indexQjob(indexJob)=ii indexForce(indexJob)=force indexTemperature(indexJob)=temperature indexNeck(indexJob)=i_neck IZ_JOB(indexJob)=zNumber IN_JOB(indexJob)=nNumber If(x_proton_initValue.Gt.0.0D0) Then XIZJOB(indexJob)=xzNumber Else XIZJOB(indexJob)=Real(zNumber,Kind=Kind(1.0)) End If If(x_neutron_initValue.Gt.0.0D0) Then XINJOB(indexJob)=xnNumber Else XINJOB(indexJob)=Real(nNumber,Kind=Kind(1.0)) End If idxDefo=mapping_qlm(ii) idxNeck=mapping_neck(i_neck) iGrid=mapping_oldGrid(idxDefo, idxNeck) indexRestartJob(indexJob)=indexJobVector_orig(iGrid, iCommon) End Do End Do End Do End Do End Do End Do Deallocate(indexJobVector_orig,mapping_oldGrid,mapping_qlm,mapping_neck) End If End If If(modeMPI_deformation.Ge.3) Then NOFJOB = proton_stepNumber * neutron_stepNumber * counterForce * numberQ & * counterTemperature * neck_numberPoints Allocate(IZ_JOB(1:NOFJOB),IN_JOB(1:NOFJOB),indexQjob(1:NOFJOB),indexRestartJob(1:NOFJOB),& indexForce(1:NOFJOB),indexTemperature(1:NOFJOB),indexNeck(1:NOFJOB)) Allocate(XIZJOB(1:NOFJOB),XINJOB(1:NOFJOB)); XIZJOB(:)=-1.0D0; XINJOB(:)=-1.0D0 nCommon = proton_stepNumber * neutron_stepNumber * counterForce * counterTemperature iCommon = 0; indexJob = 0 Do iz=1,proton_stepNumber zNumber = zVector(iz) If(x_proton_initValue.Gt.0.0D0) xzNumber=xzVector(iz) Do in=1,neutron_stepNumber nNumber = nVector(in) If(x_neutron_initValue.Gt.0.0D0) xnNumber=xnVector(in) Do force=1,counterForce iCommon = iCommon + 1 Do temperature=1,counterTemperature Do i_neck=1,neck_numberPoints Do ii=1,numberQ indexJob=indexJob+1 indexQjob(indexJob)=ii indexForce(indexJob)=force indexTemperature(indexJob)=temperature indexNeck(indexJob)=i_neck IZ_JOB(indexJob)=zNumber IN_JOB(indexJob)=nNumber If(x_proton_initValue.Gt.0.0D0) Then XIZJOB(indexJob)=xzNumber Else XIZJOB(indexJob)=Real(zNumber,Kind=Kind(1.0)) End If If(x_neutron_initValue.Gt.0.0D0) Then XINJOB(indexJob)=xnNumber Else XINJOB(indexJob)=Real(nNumber,Kind=Kind(1.0)) End If indexRestartJob(indexJob)=mapping_qlm_new(ii) If(debug.Ge.4 .And. mpi_rank .Eq.0) Then Do iq=1,numberConstraints Write(6,'("Z=",i3," N=",i3," force=",a4," T=",f8.4, & & " iq=",i2," l=",i2, & & " m=",i2," Qlm=",f10.5," -- RESTART i=",i6)') & zNumber,nNumber, & forceVector(indexForce(indexJob)), & temperatureVector(temperature), & iq,lambdaJob(indexQjob(indexJob),iq), & miuJob(indexQjob(indexJob),iq), & qValueJob(indexQjob(indexJob),iq), & ii End Do End If End Do End Do End Do End Do End Do End Do Deallocate(mapping_qlm_new,mapping_neck) End If If(modeMPI_basis==1) Then Allocate(basis_size(1:basis_size_stepNumber)) Allocate(basis_states(1:basis_states_stepNumber)) Allocate(basis_freq(1:basis_frequency_stepNumber)) Allocate(basis_defo(1:basis_deform_stepNumber)) Do i_shel=1,basis_size_stepNumber basis_size(i_shel) = basis_size_initValue + (i_shel-1)*basis_size_stepValue End do Do i_stat=1,basis_states_stepNumber basis_states(i_stat) = basis_states_initValue + (i_stat-1)*basis_states_stepValue End do do i_freq=1,basis_frequency_stepNumber basis_freq(i_freq) = basis_frequency_initValue + Real(i_freq-1)*basis_frequency_stepValue End do do i_defo=1,basis_deform_stepNumber basis_defo(i_defo) = basis_deform_initValue + Real(i_defo-1)*basis_deform_stepValue End do NOFJOB = proton_stepNumber * neutron_stepNumber * numberQ & * basis_frequency_stepNumber * basis_deform_stepNumber & * basis_size_stepNumber * basis_states_stepNumber * counterForce Allocate(IZ_JOB(1:NOFJOB),IN_JOB(1:NOFJOB),indexQjob(1:NOFJOB),indexForce(1:NOFJOB),& indexSizeJob(1:NOFJOB),indexStatesJob(1:NOFJOB),indexLengthJob(1:NOFJOB), & idxDeformJob(1:NOFJOB)) Allocate(XIZJOB(1:NOFJOB),XINJOB(1:NOFJOB)); XIZJOB(:)=-1.0D0; XINJOB(:)=-1.0D0 indexJob = 0 Do iz=1,proton_stepNumber zNumber = zVector(iz) If(x_proton_initValue.Gt.0.0D0) xzNumber=xzVector(iz) Do in=1,neutron_stepNumber nNumber = nVector(in) If(x_neutron_initValue.Gt.0.0D0) xnNumber=xnVector(in) Do force=1,counterForce Do ii=1,numberQ Do i_shel=1,basis_size_stepNumber Do i_stat=1,basis_states_stepNumber Do i_freq=1,basis_frequency_stepNumber Do i_defo=1,basis_deform_stepNumber indexJob = indexJob + 1 IZ_JOB(indexJob)=zNumber IN_JOB(indexJob)=nNumber If(x_proton_initValue.Gt.0.0D0) Then XIZJOB(indexJob)=xzNumber Else XIZJOB(indexJob)=Real(zNumber,Kind=Kind(1.0)) End If If(x_neutron_initValue.Gt.0.0D0) Then XINJOB(indexJob)=xnNumber Else XINJOB(indexJob)=Real(nNumber,Kind=Kind(1.0)) End If indexQjob(indexJob)=ii indexForce(indexJob)=force indexSizeJob(indexJob)=i_shel indexStatesJob(indexJob)=i_stat indexLengthJob(indexJob)=i_freq idxDeformJob(indexJob)=i_defo End Do End Do End Do End Do End Do End Do End Do End Do End If End Subroutine mpi_constructJobList !---------------------------------------------------------------------! ! Subroutine setting local (process-wise) HFODD variables based on ! ! existing list of tasks. Routine executed by each process. ! !---------------------------------------------------------------------! Subroutine mpi_getCurrentData(INDJOB,IN_FIX,IZ_FIX,SKYRME,GOGNAM, & INSIQN,IPSIQN,IDSIQN, & INSIQP,IPSIQP,IDSIQP, & INSIGN,IPSIGN,ISSIGN,IDSIGN, & INSIGP,IPSIGP,ISSIGP,IDSIGP, & ICONTI,IPCONT,ILCONT,IFCONT, & I_GOGA, & SLOWEV,SLOWOD,SLOWPA,SLOWLI, & IROTAT, & NILXYZ, & N_EXIT, & NOITER, & FCHOM0,NOSCIL,NLIMIT, & LIPKIN,LIPKIP, & NXHERM,NYHERM,NZHERM,IOPTGS, & IPAIRI,IPAHFB,IMFHFB, & IWRILI, & EPSITE, & INIROT, & IGAMMA, & LANODD, & IFIBLN,INIBLN,IFIBLP,INIBLP, & NMUCON, & IFSHEL, & IF_RPA, & INDFIL, & TEMP_T) Use hfodd_sizes Integer :: INDJOB,IN_FIX,IZ_FIX,INSIQN,IPSIQN,IDSIQN,INSIQP,IPSIQP,IDSIQP,& INSIGN,IPSIGN,ISSIGN,IDSIGN,INSIGP,IPSIGP,ISSIGP,IDSIGP,ICONTI,& IPCONT,ILCONT,IROTAT,NILXYZ,N_EXIT,NOITER,NOSCIL,NLIMIT,LIPKIN,& LIPKIP,NXHERM,NYHERM,NZHERM,IOPTGS,IPAIRI,IPAHFB,IMFHFB,IWRILI,& IFCONT,INIROT,IGAMMA,LANODD,IFIBLN,INIBLN,IFIBLP,INIBLP,NMUCON,& IFSHEL,IF_RPA,INDFIL,I_GOGA Integer :: IFLAGQ,IFLALQ,indexQ,ii,LAMBDA,MIU,IFNECK Integer :: i_shel,i_stat,i_freq,i_defo,force,temperature,i_neck,is_deformed Integer :: mpi_rank,mpi_err Character(Len=4) :: SKYRME,GOGNAM Real :: SLOWEV,SLOWOD,SLOWPA,SLOWLI,EPSITE,FCHOM0,STIFFQ, & QASKED,GALMUQ,QLINEA,Q0NECK,G_NECK,HBMASS,HBMRPA, & HBMINP,H_BARC,XMASSP,XMASSN,HBCOE2,HBAROX,HBAROY, & HBAROZ,HOMSCA,ALPHAR,TEMP_T,PIARGU,OVALUE,A_MASS, & HOMEGA,XINFIX,XIZFIX COMMON & /FLOPAR/ XIZFIX,XINFIX COMMON & /NCKVAL/ Q0NECK,G_NECK & /NCKFLA/ IFNECK COMMON & /REALPH/ ALPHAR(0:NDLAMB,0:NDLAMB) COMMON & /QCNSTR/ STIFFQ(0:NDMULT,-NDMULT:NDMULT), & QASKED(0:NDMULT,-NDMULT:NDMULT), & IFLAGQ(0:NDMULT,-NDMULT:NDMULT) COMMON & /QLASTR/ GALMUQ(0:NDMULT,-NDMULT:NDMULT), & QLINEA(0:NDMULT,-NDMULT:NDMULT), & IFLALQ(0:NDMULT,-NDMULT:NDMULT) COMMON & /PLANCK/ HBMASS,HBMRPA,HBMINP COMMON & /PHYCON/ H_BARC,XMASSP,XMASSN,HBCOE2 COMMON & /BASPAR/ HBAROX,HBAROY,HBAROZ & /SCALNG/ HOMSCA(NDKART) Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err) ! Set proton, neutron number and Skyrme functional IN_FIX=IN_JOB(INDJOB); IZ_FIX=IZ_JOB(INDJOB) XINFIX=XINJOB(INDJOB); XIZFIX=XIZJOB(INDJOB) If(debug.Ge.1) Write(6,'("In mpi_getCurrentData - mpi_rank = ",i10," XN = ",f20.14," XZ = ",f20.14)') & mpi_rank, XINFIX, XIZFIX If(I_GOGA.Gt.0) Then force=indexForce(INDJOB) GOGNAM=forceVector(force) SKYRME=GOGNAM Else force=indexForce(INDJOB) SKYRME=forceVector(force) End If ! Set temperature and index of the restart file (assumed to be in restart/ directory, ! and to have the generic form HFODD_XXXXXXXX.REC, where XXXXXXXX is a 8 digit integer If(modeMPI_deformation.Ge.1.Or.modeMPI_neck.Ge.1) Then temperature=indexTemperature(INDJOB) TEMP_T=temperatureVector(temperature) INDFIL=indexRestartJob(INDJOB) If(debug.Ge.1) Write(6,'("In mpi_getCurrentData - mpi_rank = ",i10," INDFIL = ",i10)') & mpi_rank,INDFIL End If ! Set constraint on the neck, and active it If(modeMPI_neck.Ge.1) Then IFNECK=2 i_neck=indexNeck(INDJOB) Q0NECK=neck_value(i_neck) End If ! Set constraints on multipole moments, and activate constraints indexQ=indexQjob(INDJOB) If(debug.Ge.1) Write(6,'("In mpi_getCurrentData - mpi_rank = ",i10," INDJOB = ",i10," indexQ = ",i10)') & mpi_rank,INDJOB,indexQ Do ii=1,numberConstraints LAMBDA=Abs(lambdaJob(indexQ,ii)) If(LAMBDA.Ge.NMUCON) NMUCON=LAMBDA MIU = miuJob(indexQ, ii) QASKED(LAMBDA,MIU) = qValueJob(indexQ, ii) IFLAGQ(LAMBDA,MIU) = 0 If(turnConstraintOn.Eq.1) IFLAGQ(LAMBDA,MIU) = 1 IFLALQ(LAMBDA,MIU) = 0 If(IF_RPA.Eq.1) IFLALQ(LAMBDA,MIU) =-1 ! RPA STIFFQ(LAMBDA,MIU) = 0.5D0 End Do ! Set characteristics of the basis for this process If(modeMPI_basis.Eq.1) Then i_shel=indexSizeJob(INDJOB) i_stat=indexStatesJob(INDJOB) i_freq=indexLengthJob(INDJOB) i_defo=idxDeformJob(INDJOB) ! Number of shells and number of states NOSCIL = basis_size(i_shel) NLIMIT = basis_states(i_stat) ! Check if the requested configuration is deformed is_deformed = 0 Do lambda=0,NDMULT Do miu=0,NDMULT If(QASKED(LAMBDA,MIU).Gt.1.D-12) Then If(debug.Ge.4) Write(6,'("LAMBDA = ",i2," MIU = ",i2," QASKED = ",f20.14)') & LAMBDA,MIU,QASKED(LAMBDA,MIU) is_deformed = 1 End If End Do End Do ! Set frequency and deformation of the basis A_MASS=Real(IZ_FIX+IN_FIX); HOMEGA=41.0D0/A_MASS**(1.0D0/3.0D0) OVALUE=basis_freq(i_freq); FCHOM0=OVALUE/HOMEGA ALPHAR(2,0)=basis_defo(i_defo) If(debug.Ge.1) Write(6,'("NOSCIL=",i4," NLIMIT=",i4," FCHOM0=",f10.5," ALPHAR=",f10.5)') & NOSCIL,NLIMIT,FCHOM0,ALPHAR(2,0) End If End Subroutine mpi_getCurrentData !---------------------------------------------------------------------! ! Subroutine that allows to redefine the requested values of the ! ! constraints based on the values of Qlm read from disk. This is ! ! used in restart mode only to effectively constrain previously free ! ! collective variables. ! ! Example: Run 1 - Q20 grid, Q40 unconstrained ! ! Run 2 - Same Q20 grid, new Q40 grid between 0 and 50. ! ! Calling mpi_pathExploration will recenter the Q40 ! ! grid around the equilibrium values of Q40 from ! ! Run 1. ! !---------------------------------------------------------------------! Subroutine mpi_pathExploration(IERROR,NMUCON) Use hfodd_sizes Integer :: IERROR,NMUCON,ITITLE,NFIPRI,IFLAGQ,LAMACT,MIUACT,i,lambda,miu Real :: QOLDIE,QSHIFT,STIFFQ,QASKED,QMUL_T,QMUL_P,QMUL_N COMMON & /CFIPRI/ NFIPRI COMMON & /QCNSTR/ STIFFQ(0:NDMULT,-NDMULT:NDMULT), & QASKED(0:NDMULT,-NDMULT:NDMULT), & IFLAGQ(0:NDMULT,-NDMULT:NDMULT) COMMON & /QMULTI/ QMUL_N(0:NDMULT,-NDMULT:NDMULT), & QMUL_P(0:NDMULT,-NDMULT:NDMULT), & QMUL_T(0:NDMULT,-NDMULT:NDMULT) If(modeMPI_deformation.Eq.2.And.optimizePath.Eq.1.And.IERROR.Eq.0) Then ! Readjustement only for extra constraints not included in the original run Do i=1,numberConstraints LAMACT=qLambda_restart(i) MIUACT=qMiu_restart(i) QOLDIE=QASKED(LAMACT,MIUACT) ! Original value requested QSHIFT=QMUL_T(LAMACT,MIUACT) - qinitValue_restart(i) & -0.5D0*(qFinValue_restart(i) - qinitValue_restart(i)) QASKED(LAMACT,MIUACT)=QSHIFT+QASKED(LAMACT,MIUACT) ! new value centered on unconstrained result End Do ! Displaying a message with the new values of the constraints Write(NFIPRI,'("*",77X,"*")') ITITLE=0 Do lambda=1,NMUCON Do miu=-lambda,lambda If(IFLAGQ(lambda,miu).Eq.1) Then If(ITITLE.Eq.0) Then Write(NFIPRI,'(79("*"),/,"*",77X,"*")') Write(NFIPRI,'("*",1X,"UPDATED MULTIPOLE CONSTRAINTS: ",& & "LAMBDA=",I2," MIU=",I2,2X,"MOMENT=",F7.3,12X,"*")') & lambda,miu,QASKED(lambda,miu) ITITLE=1 Else Write(NFIPRI,'("*",32X,"LAMBDA=",I2," MIU=",I2,2X,"MOMENT=",F7.3,12X,"*")') & lambda,miu,QASKED(lambda,miu) End If End If End Do End Do End If End Subroutine mpi_pathExploration !---------------------------------------------------------------------! ! Subroutine reading the HFODD record file. It is a modified version ! ! of routine RECORD(), where storage of big arrays is made through ! ! allocatable Fortran arrays. ! !---------------------------------------------------------------------! Subroutine read_hfodd(NFIREP,FILREP,NUMITE,NXHERM,NYHERM,NZHERM,IPCONT, & ILCONT,IACONT,IMCONT,IRCONT,IRENMA,IRENIN,LIPKIN, & LIPKIP,REFERN,REFERP,REDELN,REDELP,REFE2N,REFE2P, & IFIBLN,INIBLN,IFIBLP,INIBLP,NMUMAX,NSIMAX,NMUCON, & ISHIFT,IERROR) Use hfodd_sizes Integer :: currentVersion, minimalVersion, mpi_rank, mpi_err, iq Integer :: IX,IY,IZ,K,L,LDBASE,NFIPRI,IFNECK,IFIRED,LAMBDA,MIU, & NFIREP,NUMITE,NXHERM,NYHERM,NZHERM,MDMULT,IPCONT, & ILCONT,IACONT,IMCONT,IRCONT,IRENMA,IRENIN,LIPKIN, & LIPKIP,IFIBLN,INIBLN,IFIBLP,INIBLP,NMUMAX,NSIMAX, & NMUCON,ISHIFT,IERROR,IVERIN,MXHERM,MYHERM,MZHERM, & MMUMAX,MSIMAX,MDBASE,ICHARG,IREVER,KARTEZ,IBASE,ii Integer, Allocatable :: IFLALQ(:,:),IFLALS(:,:),IDEBLO(:) character(Len=68) :: FILREP Real :: REFERN,REFERP,REDELN,REDELP,REFE2N,REFE2P,Q0NECK,G_NECK,DUMMY Real, Allocatable :: ANGU_N(:),ANGU_P(:),ANGU_T(:),SPIN_N(:),SPIN_P(:),SPIN_T(:) Real, Allocatable :: QMUL_N(:,:),QMUL_P(:,:),QMUL_T(:,:),SMUL_N(:,:),SMUL_P(:,:), & SMUL_T(:,:),GALMUQ(:,:),QLINEA(:,:),GALMUS(:,:),SLINEA(:,:), & RALMUQ(:,:),RALMUS(:,:) Real, Allocatable :: FERINI(:),DELINI(:),FE2INI(:),RSHIFT(:),HBMREN(:),ROTREN(:), & HBMREA(:),ROTREA(:) Real, Allocatable :: DLINSN(:),DLINSP(:),DLINST(:),ELINSN(:),ELINSP(:),ELINST(:), & TLINSN(:),TLINSP(:),TLINST(:),ALINLN(:),ALINLP(:),ALINLT(:), & DROTSN(:),DROTSP(:),DROTST(:),EROTSN(:),EROTSP(:),EROTST(:), & TROTSN(:),TROTSP(:),TROTST(:),AROTLN(:),AROTLP(:),AROTLT(:) Real, Allocatable :: VN_MAS(:,:,:),VN_CEN(:,:,:),VP_MAS(:,:,:),VP_CEN(:,:,:), & DPRRHO(:,:,:),DNRRHO(:,:,:) Real, Allocatable :: VN_KIS(:,:,:,:),VN_SPI(:,:,:,:),VN_CUR(:,:,:,:), & VP_KIS(:,:,:,:),VP_SPI(:,:,:,:),VP_CUR(:,:,:,:) Real, Allocatable :: VN_SOR(:,:,:,:,:),VP_SOR(:,:,:,:,:) Complex :: C_ZERO Complex, Allocatable :: DN_RHO(:,:,:),DN_TAU(:,:,:),DN_LPR(:,:,:),DN_DIV(:,:,:), & DP_RHO(:,:,:),DP_TAU(:,:,:),DP_LPR(:,:,:),DP_DIV(:,:,:),& WN_CEN(:,:,:),WP_CEN(:,:,:),WAVBLO(:,:,:) Complex, Allocatable :: DN_SPI(:,:,:,:),DN_KIS(:,:,:,:),DN_GRR(:,:,:,:),DN_LPS(:,:,:,:),& DN_ROS(:,:,:,:),DN_ROC(:,:,:,:),DN_CUR(:,:,:,:), & DP_SPI(:,:,:,:),DP_KIS(:,:,:,:),DP_GRR(:,:,:,:),DP_LPS(:,:,:,:),& DP_ROS(:,:,:,:),DP_ROC(:,:,:,:),DP_CUR(:,:,:,:) Complex, Allocatable :: DN_SCU(:,:,:,:,:),DN_DES(:,:,:,:,:),DP_SCU(:,:,:,:,:),DP_DES(:,:,:,:,:) COMMON & /DIMENS/ LDBASE COMMON & /CFIPRI/ NFIPRI COMMON & /NCKFLA/ IFNECK C_ZERO=Cmplx(0.0D0,0.0D0) currentVersion = 16; minimalVersion = 4 Call mpi_comm_rank(MPI_COMM_WORLD,mpi_rank,mpi_err) Open(UNIT=NFIREP,FILE=FILREP,STATUS='OLD',FORM='UNFORMATTED',IOSTAT=IERROR) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in opening file ",i5)') NFIREP Return End If IFIRED=0 Read(NFIREP,IOSTAT=IERROR) IVERIN If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - error in reading IVERIN")') Return End If If(.Not.Allocated(WN_CEN)) & Allocate(WN_CEN(NDXHRM,NDYHRM,NDZHRM),VN_MAS(NDXHRM,NDYHRM,NDZHRM),VN_CEN(NDXHRM,NDYHRM,NDZHRM),& WP_CEN(NDXHRM,NDYHRM,NDZHRM),VP_MAS(NDXHRM,NDYHRM,NDZHRM),VP_CEN(NDXHRM,NDYHRM,NDZHRM)) If(.Not.Allocated(VN_KIS)) & Allocate(VN_KIS(NDXHRM,NDYHRM,NDZHRM,NDKART),VN_SPI(NDXHRM,NDYHRM,NDZHRM,NDKART),& VN_CUR(NDXHRM,NDYHRM,NDZHRM,NDKART),VP_KIS(NDXHRM,NDYHRM,NDZHRM,NDKART),& VP_SPI(NDXHRM,NDYHRM,NDZHRM,NDKART),VP_CUR(NDXHRM,NDYHRM,NDZHRM,NDKART)) If(.Not.Allocated(VN_SOR)) & Allocate(VN_SOR(NDXHRM,NDYHRM,NDZHRM,NDKART,NDKART),VP_SOR(NDXHRM,NDYHRM,NDZHRM,NDKART,NDKART)) If(.Not.Allocated(DN_RHO)) & Allocate(DN_RHO(NDXHRM,NDYHRM,NDZHRM),DN_TAU(NDXHRM,NDYHRM,NDZHRM),DN_LPR(NDXHRM,NDYHRM,NDZHRM),& DN_DIV(NDXHRM,NDYHRM,NDZHRM),DP_RHO(NDXHRM,NDYHRM,NDZHRM),DP_TAU(NDXHRM,NDYHRM,NDZHRM),& DP_LPR(NDXHRM,NDYHRM,NDZHRM),DP_DIV(NDXHRM,NDYHRM,NDZHRM)) If(.Not.Allocated(DN_SPI)) & Allocate(DN_SPI(NDXHRM,NDYHRM,NDZHRM,NDKART),DN_KIS(NDXHRM,NDYHRM,NDZHRM,NDKART),& DN_GRR(NDXHRM,NDYHRM,NDZHRM,NDKART),DN_LPS(NDXHRM,NDYHRM,NDZHRM,NDKART),& DN_ROS(NDXHRM,NDYHRM,NDZHRM,NDKART),DN_ROC(NDXHRM,NDYHRM,NDZHRM,NDKART),& DN_CUR(NDXHRM,NDYHRM,NDZHRM,NDKART)) If(.Not.Allocated(DP_SPI)) & Allocate(DP_SPI(NDXHRM,NDYHRM,NDZHRM,NDKART),DP_KIS(NDXHRM,NDYHRM,NDZHRM,NDKART),& DP_GRR(NDXHRM,NDYHRM,NDZHRM,NDKART),DP_LPS(NDXHRM,NDYHRM,NDZHRM,NDKART),& DP_ROS(NDXHRM,NDYHRM,NDZHRM,NDKART),DP_ROC(NDXHRM,NDYHRM,NDZHRM,NDKART),& DP_CUR(NDXHRM,NDYHRM,NDZHRM,NDKART)) If(.Not.Allocated(DN_SCU)) & Allocate(DN_SCU(NDXHRM,NDYHRM,NDZHRM,NDKART,NDKART),DN_DES(NDXHRM,NDYHRM,NDZHRM,NDKART,NDKART),& DP_SCU(NDXHRM,NDYHRM,NDZHRM,NDKART,NDKART),DP_DES(NDXHRM,NDYHRM,NDZHRM,NDKART,NDKART)) If(.Not.Allocated(ANGU_N)) & Allocate(ANGU_N(1:NDKART),ANGU_P(1:NDKART),ANGU_T(1:NDKART),& SPIN_N(1:NDKART),SPIN_P(1:NDKART),SPIN_T(1:NDKART)) If(.Not.Allocated(QMUL_N)) & Allocate(QMUL_N(0:NDMULT,-NDMULT:NDMULT),QMUL_P(0:NDMULT,-NDMULT:NDMULT),QMUL_T(0:NDMULT,-NDMULT:NDMULT),& SMUL_N(0:NDMULT,-NDMULT:NDMULT),SMUL_P(0:NDMULT,-NDMULT:NDMULT),SMUL_T(0:NDMULT,-NDMULT:NDMULT),& GALMUQ(0:NDMULT,-NDMULT:NDMULT),QLINEA(0:NDMULT,-NDMULT:NDMULT),RALMUQ(0:NDMULT,-NDMULT:NDMULT),& GALMUS(0:NDMULT,-NDMULT:NDMULT),SLINEA(0:NDMULT,-NDMULT:NDMULT),RALMUS(0:NDMULT,-NDMULT:NDMULT)) If(.Not.Allocated(IFLALQ)) & Allocate(IFLALQ(0:NDMULT,-NDMULT:NDMULT),IFLALS(0:NDMULT,-NDMULT:NDMULT)) If(.Not.Allocated(FERINI)) & Allocate(FERINI(0:NDISOS),DELINI(0:NDISOS),FE2INI(0:NDISOS),RSHIFT(1:NDKART),& HBMREN(1:NDKART),ROTREN(1:NDKART),HBMREA(1:NDKART),ROTREA(1:NDKART)) If(.Not.Allocated(DLINSN)) & Allocate(DLINSN(0:NDKART),DLINSP(0:NDKART),DLINST(0:NDKART), & ELINSN(0:NDKART),ELINSP(0:NDKART),ELINST(0:NDKART), & TLINSN(0:NDKART),TLINSP(0:NDKART),TLINST(0:NDKART), & ALINLN(0:NDKART),ALINLP(0:NDKART),ALINLT(0:NDKART), & DROTSN(0:NDKART),DROTSP(0:NDKART),DROTST(0:NDKART), & EROTSN(0:NDKART),EROTSP(0:NDKART),EROTST(0:NDKART), & TROTSN(0:NDKART),TROTSP(0:NDKART),TROTST(0:NDKART), & AROTLN(0:NDKART),AROTLP(0:NDKART),AROTLT(0:NDKART)) If(.Not.Allocated(WAVBLO)) & Allocate(WAVBLO(1:NDBASE,0:NDISOS,0:NDREVE)) If(.Not.Allocated(IDEBLO)) & Allocate(IDEBLO(0:NDISOS)) If(IVERIN.Ge.minimalVersion.AND.IVERIN.Le.currentVersion) Then Read(NFIREP,IOSTAT=IERROR) NUMITE If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading NUMITE")') Return End If Read(NFIREP,IOSTAT=IERROR) MXHERM,MYHERM,MZHERM If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading MXHERM")') Return End If If(MXHERM.Ne.NXHERM.OR.MYHERM.Ne.NYHERM.OR.MZHERM.Ne.NZHERM) Then Write(NFIPRI,'(1X,19("/")," WRONG NXHERM,NYHERM,NZHERM:",3I3,2X,19("/"))') MXHERM,MYHERM,MZHERM Stop ' WRONG NXHERM,NYHERM,NZHERM IN RECORD' End If Read(NFIREP,IOSTAT=IERROR) (((VN_MAS(IX,IY,IZ),VN_CEN(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading VN_MAS")') Return End If Read(NFIREP,IOSTAT=IERROR) ((((VN_KIS(IX,IY,IZ,L),VN_SPI(IX,IY,IZ,L),VN_CUR(IX,IY,IZ,L), & IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM),L=1,NDKART) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading VN_KIS")') Return End If Read(NFIREP,IOSTAT=IERROR) (((((VN_SOR(IX,IY,IZ,L,K),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM),L=1,NDKART),& K=1,NDKART) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading VN_SOR")') Return End If Read(NFIREP,IOSTAT=IERROR) (((VP_MAS(IX,IY,IZ),VP_CEN(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading VP_MAS")') Return End If Read(NFIREP,IOSTAT=IERROR) ((((VP_KIS(IX,IY,IZ,L),VP_SPI(IX,IY,IZ,L),VP_CUR(IX,IY,IZ,L), & IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM),L=1,NDKART) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading VP_KIS")') Return End If Read(NFIREP,IOSTAT=IERROR) (((((VP_SOR(IX,IY,IZ,L,K),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM),L=1,NDKART),& K=1,NDKART) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading VP_SOR")') Return End If If(IVERIN.Le.6) Then If(IVERIN.Gt.5) Then Read(NFIREP,IOSTAT=IERROR) (((DNRRHO(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading DNRRHO")') Return End If Else DNRRHO(:,:,:)=0.0D0 End If Read(NFIREP,IOSTAT=IERROR) (((DPRRHO(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading DPRRHO")') Return End If DN_RHO(:,:,:)=DNRRHO(:,:,:)*cmplx(1.0d0,0.0d0) DP_RHO(:,:,:)=DPRRHO(:,:,:)*cmplx(1.0d0,0.0d0) Else Read(NFIREP,IOSTAT=IERROR) (((DN_RHO(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading DN_RHO")') Return End If Read(NFIREP,IOSTAT=IERROR)(((DP_RHO(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading DP_RHO")') Return End If End If Read(NFIREP,IOSTAT=IERROR) MMUMAX If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading NMUMAX")') Return End If If(MMUMAX.Lt.NMUMAX) Then Write(NFIPRI,'(1X,20("/")," WARNING: TOO FEW MULTIPOLE MOMENTS",1X,20("/"))') Write(NFIPRI,'(1X,20("/")," FOUND ON THE REPLAY FILE. ",1X,20("/"))') Write(NFIPRI,'(1X,20("/")," REQUESTED MAXIMUM MULTIPOLARITY=",I2,1X,20("/"))') NMUMAX Write(NFIPRI,'(1X,20("/")," FOUND MAXIMUM MULTIPOLARITY=",I2,1X,20("/"))') MMUMAX Write(NFIPRI,'(1X,20("/")," A CONSTRAINT ON A MISSING MULTIPOLE",1X,20("/"))') Write(NFIPRI,'(1X,20("/")," WILL BE INCLUDED ONLY AFTER THE",1X,20("/"))') Write(NFIPRI,'(1X,20("/")," FIRST DIAGONALIZATION IS PERFORMED",1X,20("/"))') Do LAMBDA=MMUMAX+1,NMUMAX Do MIU=-LAMBDA,LAMBDA QMUL_N(LAMBDA,MIU)=0.0D0 QMUL_P(LAMBDA,MIU)=0.0D0 QMUL_T(LAMBDA,MIU)=0.0D0 End Do End Do End If Read(NFIREP,IOSTAT=IERROR) ((QMUL_N(LAMBDA,MIU),QMUL_P(LAMBDA,MIU),QMUL_T(LAMBDA,MIU), & MIU=-LAMBDA,LAMBDA),LAMBDA=0,MMUMAX) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading QMUL_N")') Return End If Do LAMBDA=0,NSIMAX Do MIU=-LAMBDA,LAMBDA SMUL_N(LAMBDA,MIU)=0.0D0 SMUL_P(LAMBDA,MIU)=0.0D0 SMUL_T(LAMBDA,MIU)=0.0D0 End Do End Do If(IVERIN.Gt.4) Then Read(NFIREP,IOSTAT=IERROR) MSIMAX If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading MSIMAX")') Return End If If(MSIMAX.Lt.NSIMAX) Then Write(NFIPRI,'(1X,20(1H/)," WARNING: TOO FEW SURFACE MOMENTS",1X,20(1H/))') Write(NFIPRI,'(1X,20(1H/)," FOUND ON THE REPLAY FILE. ",1X,20(1H/))') Write(NFIPRI,'(1X,20(1H/)," REQUESTED MAXIMUM MULTIPOLARITY=",I2,1X,20(1H/))') NSIMAX Write(NFIPRI,'(1X,20(1H/)," FOUND MAXIMUM MULTIPOLARITY=",I2,1X,20(1H/))') MSIMAX Write(NFIPRI,'(1X,20(1H/)," A CONSTRAINT ON A MISSING MULTIPOLE",1X,20(1H/))') Write(NFIPRI,'(1X,20(1H/)," WILL BE INCLUDED ONLY AFTER THE",1X,20(1H/))') Write(NFIPRI,'(1X,20(1H/)," FIRST DIAGONALIZATION IS PERFORMED",1X,20(1H/))') Do LAMBDA=MSIMAX+1,NSIMAX Do MIU=-LAMBDA,LAMBDA SMUL_N(LAMBDA,MIU)=0.0D0 SMUL_P(LAMBDA,MIU)=0.0D0 SMUL_T(LAMBDA,MIU)=0.0D0 End Do End Do End If Read(NFIREP,IOSTAT=IERROR) ((SMUL_N(LAMBDA,MIU),SMUL_P(LAMBDA,MIU),SMUL_T(LAMBDA,MIU), & MIU=-LAMBDA,LAMBDA),LAMBDA=0,MSIMAX) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading SMUL_N")') Return End If End If Read(NFIREP,IOSTAT=IERROR) ISHIFT,(RSHIFT(L),L=1,NDKART) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading ISHIFT")') Return End If Read(NFIREP,IOSTAT=IERROR) (ANGU_N(L),ANGU_P(L),ANGU_T(L),SPIN_N(L),SPIN_P(L),SPIN_T(L),L=1,NDKART) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading ANGU_N")') Return End If Read(NFIREP,IOSTAT=IERROR) (((WN_CEN(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading WN_CEN")') Return End If Read(NFIREP,IOSTAT=IERROR) (((WP_CEN(IX,IY,IZ),IX=1,NXHERM),IY=1,NYHERM),IZ=1,NZHERM) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading WP_CEN")') Return End If Read(NFIREP,IOSTAT=IERROR) REFERN,REDELN,REFERP,REDELP If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading REFERN")') Return End If ! LN lambda_2 (IVERIN=8) If(IVERIN.Gt.7) Then Read(NFIREP,IOSTAT=IERROR) REFE2N,REFE2P If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading REFE2N")') Return End If Else REFE2N=FE2INI(0); REFE2P=FE2INI(1) If((LIPKIN.EQ.1.OR.LIPKIP.EQ.1).AND.ILCONT.EQ.1) Then Write(NFIPRI,'(1X,13("/")," WARNING: THE LIPKIN-NOGAMI LAMBDA2S HAVE NOT BEEN",1X,13("/"))') Write(NFIPRI,'(1X,13("/")," FOUND ON THE OLD REPLAY FILE VERSION NO. =",I2,1X,13("/"))') IVERIN Write(NFIPRI,'(1X,13("/")," THE LN CORRECTIONS CANNOT BE SMOOTHLY CONTINUED",1X,13("/"))') End If End If ! Blocked state wave-function (IVERIN=9) IDEBLO(0)=0; IDEBLO(1)=0 If(IVERIN.Gt.8) Then Read(NFIREP,IOSTAT=IERROR) MDBASE,(IDEBLO(ICHARG),ICHARG=0,NDISOS) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading MDBASE")') Return End If If(MDBASE.Le.NDBASE) Then Read(NFIREP,IOSTAT=IERROR) (((WAVBLO(IBASE,ICHARG,IREVER),IBASE=1,MDBASE),ICHARG=0,NDISOS),& IREVER=0,NDREVE) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading WAVBLO")') Return End If Else Read(NFIREP,IOSTAT=IERROR) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading ...nothing...")') Return End If End If If(((IFIBLN.EQ.1.AND.INIBLN.EQ.0).OR.(IFIBLP.EQ.1.AND.INIBLP.EQ.0)).AND. & (MDBASE.Ne.LDBASE.OR.MDBASE.Gt.NDBASE)) Then Write(NFIPRI,'(1X,10("/")," INCOMPATIBLE DIMENSIONS FOUND ON THE RECORD FILE:",1X,10("/"))') Write(NFIPRI,'(1X,10("/")," LDBASE=",I4," MDBASE=",I4," NDBASE=",I4,1X,10("/"))') LDBASE,MDBASE,NDBASE Stop 'MDBASE<>LDBASE OR MDBASE>NDBASE IN RECORD' End If End If If(IVERIN.Ge.10) Then Read(NFIREP,IOSTAT=IERROR) (ALINLT(KARTEZ),KARTEZ=0,NDKART) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading ALINLT")') Return End If Else ALINLT(:)=0.0D0 If(IMCONT.EQ.1) Then Write(NFIPRI,'(1X,13(1H/)," WARNING: THE AVERAGE LINEAR MOMENTA HAVE NOT BEEN",1X,13(1H/))') Write(NFIPRI,'(1X,13(1H/)," FOUND ON THE OLD REPLAY FILE VERSION NO. =",I2,1X,13(1H/))') IVERIN Write(NFIPRI,'(1X,13(1H/)," THE EXACT CM CORREC. CANNOT BE SMOOTHLY CONTINUED",1X,13(1H/))') End If End If If(IVERIN.Ge.11) Then Read(NFIREP,IOSTAT=IERROR) (AROTLT(KARTEZ),KARTEZ=0,NDKART) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading AROTLT")') Return End If Else AROTLT(:)=0.0D0 If(IRCONT.Eq.1) Then Write(NFIPRI,'(1X,13(1H/)," WARNING: THE AVERAGE ANGULAR MOMENTA HAVE NOT BEEN",1X,13(1H/))') Write(NFIPRI,'(1X,13(1H/)," FOUND ON THE OLD REPLAY FILE VERSION NO. =",I2,1X,13(1H/))') IVERIN Write(NFIPRI,'(1X,13(1H/)," THE EXACT ROT CORREC. CANNOT BE SMOOTHLY CONTINUED",1X,13(1H/))') End If End If If(IVERIN.Ge.12) Then Read(NFIREP,IOSTAT=IERROR) MDMULT If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading MDMULT")') Return End If RALMUQ(:,:)=0.0D0; RALMUS(:,:)=0.0D0 If(MDMULT.Ne.NDMULT) Then Write(NFIPRI,'(1X,20(1H/)," WARNING: NOT ALL LINEAR CONSTRAINTS",1X,20(1H/))') Write(NFIPRI,'(1X,20(1H/)," FOUND ON THE REPLAY FILE. ",1X,20(1H/))') Write(NFIPRI,'(1X,20(1H/)," REQUESTED MAXIMUM MULTIPOLARITY=",I2,1X,20(1H/))') NDMULT Write(NFIPRI,'(1X,20(1H/)," FOUND MAXIMUM MULTIPOLARITY=",I2,1X,20(1H/))') MDMULT Write(NFIPRI,'(1X,20(1H/)," A CONSTRAINT ON A MISSING MULTIPOLE",1X,20(1H/))') Write(NFIPRI,'(1X,20(1H/)," WILL BE INCLUDED ONLY AFTER THE",1X,20(1H/))') Write(NFIPRI,'(1X,20(1H/)," FIRST DIAGONALIZATION IS PERFORMED",1X,20(1H/))') RALMUQ(:,:)=QLINEA(:,:); RALMUS(:,:)=SLINEA(:,:) Read(NFIREP,IOSTAT=IERROR) ((DUMMY,DUMMY,MIU=-LAMBDA,LAMBDA),LAMBDA=0,MDMULT) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading DUMMY")') Return End If Else Read(NFIREP,IOSTAT=IERROR) ((RALMUQ(LAMBDA,MIU),RALMUS(LAMBDA,MIU),MIU=-LAMBDA,LAMBDA),& LAMBDA=0,MDMULT) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading RALMUQ")') Return End If ! If not allocated, allocate; if already allocated, deallocate first and re-allocate ( ! (enables vayring number of constraints) If(Allocated(flex_lambda)) Deallocate(flex_lambda,flex_miu,flex_qlm) Allocate(flex_lambda(numberConstraints),flex_miu(numberConstraints),flex_qlm(numberConstraints)) Do iq=1,numberConstraints lambda = qLambda(iq); miu = qMiu(iq) flex_lambda(iq)=lambda flex_miu(iq)=miu flex_qlm(iq)=QMUL_T(lambda,miu) End Do End If Else RALMUQ(:,:)=QLINEA(:,:); RALMUS(:,:)=SLINEA(:,:) If(IACONT.EQ.1) Then Write(NFIPRI,'(1X,13(1H/)," WARNING: THE LINEAR CONSTRAINTS H A V E NOT BEEN",1X,13(1H/))') Write(NFIPRI,'(1X,13(1H/)," FOUND ON THE OLD REPLAY FILE VERSION NO. =",I2,1X,13(1H/))') IVERIN Write(NFIPRI,'(1X,13(1H/)," THE EXACT AUGM.LAGRA. CANNOT BE SMOOTHLY CONTINUED",1X,13(1H/))') End If End If If(IVERIN.Ge.13) Then Read(NFIREP,IOSTAT=IERROR) (HBMREA(KARTEZ),KARTEZ=1,NDKART) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading HBMREA")') Return End If Read(NFIREP,IOSTAT=IERROR) (ROTREA(KARTEZ),KARTEZ=1,NDKART) If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading ROTREA")') Return End If If(IMCONT.EQ.1.AND.IRENMA.Ge.1) Then HBMREN(:)=HBMREA(:) End If If(IRCONT.EQ.1.AND.IRENIN.Ge.1) Then ROTREN(:)=ROTREA(:) End If Else If(IMCONT.EQ.1.AND.IRENMA.Ge.1) Then Write(NFIPRI,'(1X,13(1H/)," WARNING: THE RENORMALIZED MASSES HAVE NOT BEEN",1X,13(1H/))') Write(NFIPRI,'(1X,13(1H/)," FOUND ON THE OLD REPLAY FILE VERSION NO. =",I2,1X,13(1H/))') IVERIN Write(NFIPRI,'(1X,13(1H/)," THE EXACT CM CORREC. CANNOT BE SMOOTHLY CONTINUED",1X,13(1H/))') End If If(IRCONT.EQ.1.AND.IRENIN.Ge.1) Then Write(NFIPRI,'(1X,13(1H/)," WARNING: THE RENORMALIZED INERTIA HAVE NOT BEEN",1X,13(1H/))') Write(NFIPRI,'(1X,13(1H/)," FOUND ON THE OLD REPLAY FILE VERSION NO. =",I2,1X,13(1H/))') IVERIN Write(NFIPRI,'(1X,13(1H/)," THE EXACT ROT CORREC. CANNOT BE SMOOTHLY CONTINUED",1X,13(1H/))') End If End If ! Linear Constraints for the neck (IVERIN=14) If(IVERIN.Ge.14) Then Read(NFIREP,IOSTAT=IERROR) G_NECK If(IERROR.Ne.0) Then Write(NFIPRI,'("read_hfodd() - Error in reading MDMULT")') Return End If Else G_NECK=0.0D0 If(IFNECK.Ge.2) Then Write(NFIPRI,'(1X,13(1H/)," WARNING : THE NECK CONSTRAINTS H A S NOT BEEN",1X,13(1H/))') Write(NFIPRI,'(1X,13(1H/)," FOUND ON THE OLD REPLAY FILE VERSION NO. =",I2,1X,13(1H/))') IVERIN Write(NFIPRI,'(1X,13(1H/)," THE NECK CONSTRAINT CAN NOT BE SMOOTHLY CONTINUED",1X,13(1H/))') End If End If IFIRED=1 End If If(IFIRED.EQ.0) Then Write(NFIPRI,'(/,1X,10(1H/)," CANNOT read THE REPLAY FILE VERSION:",I3,2X,10(1H/),/)') IVERIN Stop ' WRONG IVERIN IN RECORD' End If If(IFIBLN.EQ.1.AND.INIBLN.EQ.0.AND.IDEBLO(0).Ne.1) Then Write(NFIPRI,'(1X,10(1H/)," NEUTRON SINGLE-PARTICLE WF NOT FOUND ON THE RECORD FILE",1X,10(1H/))') Write(NFIPRI,'(1X,10(1H/)," RERUN THE BLOCKING CALCULATION WITH INIBLN=1 ",1X,10(1H/))') Stop ' NEUTRON SINGLE-PARTICLE WF NOT FOUND IN RECORD' End If If(IFIBLP.EQ.1.AND.INIBLP.EQ.0.AND.IDEBLO(1).Ne.1) Then Write(NFIPRI,'(1X,10(1H/)," PROTON SINGLE-PARTICLE WF NOT FOUND ON THE RECORD FILE",1X,10(1H/))') Write(NFIPRI,'(1X,10(1H/)," RERUN THE BLOCKING CALCULATION WITH INIBLP=1 ",1X,10(1H/))') Stop ' PROTON SINGLE-PARTICLE WF NOT FOUND IN RECORD' End If End Subroutine read_hfodd !---------------------------------------------------------------------! ! Subroutine reads lists all the files in 'directory_name', extracts ! ! all files of the form 'HFODD_XXXXXXXX.REC' (with XXXXXXXX a number)! ! and checks that all files with the number between 1 and N are ! ! present. It then broadcasts a boolean flag indicating if all files ! ! are there, and a vector of integers giving the breakdown. ! !---------------------------------------------------------------------! Subroutine check_directory(directory_name,N) Integer, INTENT(IN) :: N Character(Len=132), INTENT(IN) :: directory_name Integer, Parameter :: ndatin=85 Integer :: i,exitstat,iocheck,ierr,mpi_rank,mpi_err,numero,count Character(Len=80) :: filein Character(Len=150) :: commande Character(Len=132) :: string mpi_rank=0; mpi_err=0 ! Getting the rank of the current process Call mpi_comm_rank(MPI_COMM_WORLD, mpi_rank, mpi_err) If(mpi_rank.Eq.0) Then commande = 'ls '//Trim(directory_name)// ' > toto'; exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("mpi_rank = ",i10," - Executed command: ",a)') & mpi_rank,commande ! open file Open(ndatin,file='toto',status='old',form='formatted',iostat=ierr) If(ierr.Ne.0) Then Write(6,'("Error in opening file : toto")') Stop 'Error in check_directory - I/O' End if ! Read file, extract number of file if relevant, and fill table iocheck=0 Do While(iocheck.Eq.0) Read(ndatin,*,Iostat=iocheck) string numero = 0 Call get_numero(string,numero) If(debug.Ge.2) Write(6,'("numero = ",i10)') numero If(numero.Gt.0) table_fichiers(numero)=1 End Do Close(ndatin) If(debug.Ge.1) Write(6,'("Done reading with the file")') filled=.True.; count=0 Do i=1,N If(table_fichiers(i).Eq.0) Then count=count+1 End If End Do If(count.Gt.0) filled=.False. If(debug.Ge.2) Then If(filled) Then Write(6,'("All files are there")') Else Write(6,'(i10," files are missing")') count End If End If End If ! Broadcast flag telling if array is filled and table of missing files Call MPI_Bcast(filled,1,MPI_LOGICAL,0,MPI_COMM_WORLD,mpi_err) Call MPI_Bcast(table_fichiers,N,MPI_INTEGER,0,MPI_COMM_WORLD,mpi_err) ! Wait for all threads Call mpi_barrier(MPI_COMM_WORLD,mpi_err) End Subroutine check_directory !---------------------------------------------------------------------! ! Subroutine returns in 'numero' the file number from a string of ! ! the form 'HFODD_XXXXXXXX.REC' (with XXXXXXXX the number). ! !---------------------------------------------------------------------! Subroutine get_numero(string,numero) Character(Len=132), INTENT(IN) :: string Integer, INTENT(INOUT) :: numero Integer :: start Character(Len=4) :: dum_4 Character(Len=6) :: dum_6 start = Index(string, 'HFODD_', .False.) If (start.Ge.0) Then Read(string(start:start+18), '(A6,I8.8,A4)') dum_6,numero,dum_4 End If End Subroutine get_numero !---------------------------------------------------------------------! ! Subroutine copies all files that have not been copied yet. ! !---------------------------------------------------------------------! Subroutine fill_holes(N,IF_THO,LIPKIN,LIPKIP) Integer, INTENT(IN) :: N,IF_THO,LIPKIN,LIPKIP Integer :: jj,exitstat,mpi_rank,mpi_world,mpi_err Character(Len=150) :: commande mpi_rank=0; mpi_world=1; mpi_err=0 ! Getting the rank of the current process Call mpi_comm_rank(MPI_COMM_WORLD, mpi_rank, mpi_err) Call mpi_comm_size(MPI_COMM_WORLD, mpi_world,mpi_err) Do jj=1,N If(Mod(mpi_rank,mpi_world).Eq.Mod(jj-1,mpi_world)) Then If(table_fichiers(jj).Eq.0) Then commande = 'cp '//fichier_rec_old(jj)//' '//fichier_rec_new(jj); exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("fill_holes: mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') & mpi_rank,jj,commande If(IF_THO.Ge.1) Then commande = 'cp '//fichier_tho_old(jj)//' '//fichier_tho_new(jj); exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("fill_holes: mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') & mpi_rank,jj,commande End If If(LIPKIN.Ge.1.Or.LIPKIP.Ge.1) Then commande = 'cp '//fichier_lic_old(jj)//' '//fichier_lic_new(jj); exitstat=0 Call system(commande) If(debug.Ge.2) Write(6,'("fill_holes: mpi_rank = ",i10," jj = ",i10," - Executed command: ",a)') & mpi_rank,jj,commande End If End If End If End Do ! Wait for all threads Call mpi_barrier(MPI_COMM_WORLD,mpi_err) End Subroutine fill_holes End Module hfodd_mpimanager