Merge pull request #117 from sharc-md/main

Merged main bugfixes into MemoryFix branch

bin/SHARC_MOLCAS.py

@@ -1542,6 +1542,11 @@ def getQMout(out, QMin):
             for jstate, j in enumerate(QMin['statemap']):
                 mult1, state1, ms1 = tuple(QMin['statemap'][i])
                 mult2, state2, ms2 = tuple(QMin['statemap'][j])
+                if 'overlap_nacs' in QMin:
+                    if mult1 not in [i[0] for i in QMin['nacmap']]:
+                        continue
+                    if mult2 not in [i[0] for i in QMin['nacmap']]:
+                        continue
                 if mult1 == mult2 and ms1 == ms2:
                     nac[istate][jstate] = complex(getsmate(out, mult1, state1, state2, states))
                 else:
@@ -2893,14 +2898,16 @@ def gettasks(QMin):
         # RASSI for overlaps
         if 'overlap' in QMin:
             if 'overlap_nacs' in QMin:
+                if imult+1 not in [ i[0] for i in QMin['nacmap']]:
+                    continue
                 tasks.append(['link', 'MOLCAS.%i.JobIph' % (imult + 1), 'JOB001'])
                 tasks.append(['link', 'MOLCAS.JobIph', 'JOB002'])
             else:
                 if 'displacement' in QMin:
                     tasks.append(['link', os.path.join(QMin['savedir'], 'MOLCAS.%i.JobIph.master' % (imult + 1)), 'JOB001'])
                 else:
                     tasks.append(['link', os.path.join(QMin['savedir'], 'MOLCAS.%i.JobIph.old' % (imult + 1)), 'JOB001'])
-            tasks.append(['link', 'MOLCAS.%i.JobIph' % (imult + 1), 'JOB002'])
+                tasks.append(['link', 'MOLCAS.%i.JobIph' % (imult + 1), 'JOB002'])
             tasks.append(['rassi', 'overlap', [nstates, nstates]])
 
         # RASSI for Dipole moments only if overlap-RASSI is not needed
@@ -3505,12 +3512,18 @@ def generate_joblist(QMin):
             QMin3['gradmap'] = [grad]
             QMin3['nacmap'] = []
             QMin3['ncpu'] = cpu_per_run[icount]
+            #for i in range(len(QMin3['states'])):
+            #    if not i+1 == grad[0]:
+            #        QMin3['states'][i] = 0
             icount += 1
             joblist[-1]['grad_%i_%i' % grad] = QMin3
         for nac in QMin['nacmap']:
             QMin3 = deepcopy(QMin2)
             QMin3['nacmap'] = [nac]
             QMin3['gradmap'] = []
+            #for i in range(len(QMin3['states'])):
+            #    if not i+1 == nac[0]:
+            #        QMin3['states'][i] = 0
             QMin3['overlap'] = [[j + 1, i + 1] for i in range(QMin['nmstates']) for j in range(i + 1)]
             QMin3['overlap_nacs'] = []
             QMin3['ncpu'] = cpu_per_run[icount]

bin/SHARC_ORCA.py

@@ -3934,14 +3934,22 @@ def get_smat_from_gbw(file1, file2=''):
             NAO = int(line.split()[0]) + 1
             break
 
+    # find start of matrix
+    iline = -1
+    while True:
+        iline += 1
+        line = out[iline]
+        if 'FRAGMENT-FRAGMENT OVERLAP MATRIX' in line:
+            break
+
     # read matrix
     nblock = 6
     ao_ovl = [[0. for i in range(NAO)] for j in range(NAO)]
     for x in range(NAO):
         for y in range(NAO):
             block = x // nblock
             xoffset = x % nblock + 1
-            yoffset = block * (NAO + 1) + y + 10
+            yoffset = block * (NAO + 1) + y + 3 + iline
             ao_ovl[x][y] = float(out[yoffset].split()[xoffset])
 
     return NAO, ao_ovl

bin/SHARC_RICC2.py

@@ -3759,12 +3759,15 @@ def define(path, QMin, ricc2=True):
 '''
     string += '''b
 all %s
-*
+''' % (QMin['template']['basis'])
+    if not ricc2:
+        string += 'c\nall 0.\n'
+    string += '''*
 eht
 y
 %i
 y
-''' % (QMin['template']['basis'],
+''' % (
        QMin['template']['charge']
        )
 

bin/diagnostics.py

@@ -755,6 +755,8 @@ def check_energies(path, trajectories, INFOS, hops):
         f = readfile(f)
         try:
             problem, tana_length = check_consistency(path, trajectories, f, 'energy.out')
+        except KeyboardInterrupt:
+            raise
         except BaseException:
             print('\n    An error occured while trying to check energy.out for consistency.')
             trajectories[path]['error'] = True
@@ -834,6 +836,8 @@ def check_populations(path, trajectories, INFOS):
         problem = ''
         try:
             problem, tana_length = check_consistency(path, trajectories, f, 'coeff_diag.out')
+        except KeyboardInterrupt:
+            raise
         except BaseException:
             print('\n    An error occured while trying to check coeff_diag.out for consistency.')
             trajectories[path]['error'] = True
@@ -901,7 +905,7 @@ def check_intruders(path, trajectories, INFOS, lis, tana, problem_length):
                 tstep = int(line.split()[3])
                 if tstep == 0:
                     prevstep = 0
-                elif tstep - 1 != prevstep:
+                elif tstep - 1 != prevstep and tstep != prevstep:
                     tana = prevstep * trajectories[path]['dtstep']
                     problem = 'Missing steps in output.log'
                     break
@@ -910,14 +914,15 @@ def check_intruders(path, trajectories, INFOS, lis, tana, problem_length):
                     ok = False
                     problem = problem_length
                     break
-            if 'RESTART requested.' in line:
-                prevstep -= 1
-            if 'time step in restart files' in line:
-                prevstep += 1
+#            if 'RESTART requested.' in line:
+#                prevstep -= 1
+#            if 'time step in restart files' in line or "STOP" in line:
+#                prevstep += 1
             if 'State: ' in line:
                 intruder = int(line.split()[1])
                 if not notpossible:
-                    state = lis[tstep][2]
+                    state = int(lis[tstep][2])
+                    # print(intruder, state)
                     if state == intruder:
                         problem = 'Intruder state found'
                         ok = False
@@ -970,6 +975,8 @@ def do_calc(INFOS):
 
             try:
                 trajectories, missing = check_files(path, trajectories, INFOS)
+            except KeyboardInterrupt:
+                raise
             except BaseException:
                 print('\n    An error occured while trying to look for the files.\n')
                 trajectories[path]['error'] = True
@@ -979,6 +986,8 @@ def do_calc(INFOS):
 
             try:
                 trajectories, f = check_runtime(path, trajectories, INFOS)
+            except KeyboardInterrupt:
+                raise
             except BaseException:
                 print('\n    An error occured while trying to extract the runtime.\n \
    Files may be corrupted.\n')
@@ -987,6 +996,8 @@ def do_calc(INFOS):
 
             try:
                 trajectories = check_termination(path, trajectories, INFOS, f)
+            except KeyboardInterrupt:
+                raise
             except BaseException:
                 print('\n    An error occured while trying to extract the status.\n \
    Files may be corrupted.\n')
@@ -1063,25 +1074,33 @@ def do_calc(INFOS):
                 step += 1
             try:
                 problem, tana = check_consistency(path, trajectories, f, 'output.lis')
+            except KeyboardInterrupt:
+               raise
             except BaseException:
                 print('\n    An error occured while trying to check output.lis for consistency.')
                 trajectories[path]['error'] = True
             if problem == '':
                 problem = check_length(path, trajectories, len(lis), 'output.lis')
             try:
                 trajectories = check_energies(path, trajectories, INFOS, hops)
+            except KeyboardInterrupt:
+                raise
             except BaseException:
                 print('\n    An error occured while trying to extract the energies.\n \
    Files may be corrupted.\n')
                 trajectories[path]['error'] = True
             try:
                 trajectories = check_populations(path, trajectories, INFOS)
+            except KeyboardInterrupt:
+                raise
             except BaseException:
                 print('\n    An error occured while trying to extract the populations.\n \
    Files may be corrupted.\n')
                 trajectories[path]['error'] = True
             try:
                 trajectories = check_intruders(path, trajectories, INFOS, lis, tana, problem)
+            except KeyboardInterrupt:
+                raise
             except BaseException:
                 print('\n    An error occured while trying to extract possible intruder states.\n \
    Files may be corrupted.\n')

bin/setup_traj.py

@@ -278,7 +278,7 @@ def all(iterable):
       'description': 'mixed gradients are calculated as linear combination of MCH gradients only',
       'required':   []
       },
-  2: {'name':        'ngh',
+  2: {'name':        'ngt',
       'description': 'mixed gradients are calculated by correction of MCH gradients with non-adiabatic coupling vector',
       'required':   ['nacdr']
      },

source/data.inc

@@ -1,2 +1,3 @@
 #define STRING_SIZE_S_ 256
 #define STRING_SIZE_L_ 1024
+#define STRING_SIZE_XL_ 4096

source/driver.f90

@@ -300,8 +300,13 @@ subroutine fixed_velocity_verlet(traj,ctrl)
     call Decoherence(traj,ctrl)
     ! obtain the correct gradient
     call Calculate_cMCH(traj,ctrl)
-    if (ctrl%calc_grad>=1) call redo_qm_gradients(traj,ctrl)
-    if (traj%kind_of_jump/=0) call Mix_gradients(traj,ctrl)
+    if (ctrl%calc_grad>=1) then
+      call redo_qm_gradients(traj,ctrl)
+      call NAC_processing(traj, ctrl)
+    endif
+    if (traj%kind_of_jump/=0) then
+      call Mix_gradients(traj,ctrl)
+    endif
 
   else if (ctrl%method==1) then !SCP
     ! Propagation coherent coefficients

source/interface.F90

@@ -464,6 +464,7 @@ subroutine get_nacdr(string, ICALL)
     use iso_c_binding
     use memory_module, only: traj, ctrl
     implicit none
+    ! __C_OUT_STRING_XL_ :: string
     __C_OUT_STRING_L_ :: string
     __INT__, intent(in) :: ICALL
     integer :: i,j
@@ -480,6 +481,7 @@ subroutine get_nacdr(string, ICALL)
             do i=1,ctrl%nstates
               do j=1,ctrl%nstates
                 if (traj%selt_ss(j,i)) write(string,'(A,1X,I3,1X,I3)') trim(string) // C_NEW_LINE , i,j
+                ! if (traj%selt_ss(j,i)) write(string,'(A,1X,I3,1X,I3)') trim(string) , i,j
               enddo
             enddo
             write(string,'(A)') trim(string) // C_NEW_LINE // 'END'
@@ -1250,7 +1252,7 @@ subroutine Verlet_finalize(IExit, iskip)
     use memory_module, only: traj, ctrl
     use misc
     use definitions
-    use qm, only: Update_old, Mix_gradients
+    use qm, only: Update_old, Mix_gradients, NAC_processing
     use electronic, only: kill_after_relaxation
     use output, only: allflush, write_dat, write_list_line, write_geom
    use restart, only: write_restart_traj
@@ -1259,7 +1261,10 @@ subroutine Verlet_finalize(IExit, iskip)
     __INT__, intent(out) :: IExit ! if IExit = 0 end loop, else continue
     __INT__, intent(in)  :: iskip ! if IExit = 0 end loop, else continue
 
-    if (traj%kind_of_jump/=0) call Mix_gradients(traj, ctrl)
+    if (traj%kind_of_jump/=0) then
+        call NAC_processing(traj, ctrl)
+        call Mix_gradients(traj, ctrl)
+    endif
     ! Finalization: Variable update, Output, Restart File, Consistency Checks
     call Update_old(traj,ctrl)
     call set_time(traj)

source/sharc_fortran.inc

@@ -5,5 +5,7 @@
 #define __COMPLEX__        complex*16
 #define __C_STRING_S_      character(len=STRING_SIZE_S_)
 #define __C_STRING_L_      character(len=STRING_SIZE_L_)
+#define __C_STRING_XL_      character(len=STRING_SIZE_XL_)
 #define __C_OUT_STRING_S_  __C_STRING_S_, intent(out) 
 #define __C_OUT_STRING_L_  __C_STRING_L_, intent(out) 
+#define __C_OUT_STRING_XL_  __C_STRING_XL_, intent(out)