From d61602b6cb4c46ba02550176b27e1ade0caae41c Mon Sep 17 00:00:00 2001
From: rpseng <rafael.enq@gmail.com>
Date: Thu, 13 Mar 2025 18:03:18 -0300
Subject: [PATCH] Conditional cosmo file generation.

---
 src/solvation/cosmo_input.F | 10 ++++++
 src/solvation/hnd_coschg.F  | 72 +++++++++++++++++++------------------
 2 files changed, 48 insertions(+), 34 deletions(-)

diff --git a/src/solvation/cosmo_input.F b/src/solvation/cosmo_input.F
index 42f272bb61..f0b9df937f 100644
--- a/src/solvation/cosmo_input.F
+++ b/src/solvation/cosmo_input.F
@@ -79,6 +79,7 @@ c
       double precision dielec,dielecinf,rsolv,rad,zeta,gammas,swtol
       character*255 field
       logical do_cosmo
+      logical do_cosmo_file
       logical do_gasphase
       logical status
       character*30 tag
@@ -326,6 +327,15 @@ c
      $                      do_cosmo_model))
      $        call errquit('cosmo_input: rtdb put failed',911, RTDB_ERR)
 c
+c   --- do_cosmo_file flag
+c
+      else if(inp_compare(.false.,'do_cosmo_file',field)) then
+        do_cosmo_file=.true.
+        if (.not. 
+     $    rtdb_put(rtdb,'cosmo:do_cosmo_file',mt_log,1,do_cosmo_file))
+     $    call errquit('cosmo_input: cannot put do_cosmo_file in rtdb',
+     $    0,rtdb_err)
+c
 c   --- rsolv: Used in Klamt-Schuurmann, but not in York-Karplus
 c   ---        approach
 c
diff --git a/src/solvation/hnd_coschg.F b/src/solvation/hnd_coschg.F
index a1bd8f0364..a98ba1ad3a 100644
--- a/src/solvation/hnd_coschg.F
+++ b/src/solvation/hnd_coschg.F
@@ -479,7 +479,7 @@ c
       double precision efcs(nefc)   !< [Input] the COSMO charge
                                     !< surface area for KS
       double precision msrf(nefc)   !< [Input] the COSMO charge
-                                    !< surface area for YS
+                                    !< surface area for YK
 c
       double precision efcz(nefc)   !< [Output] the COSMO charges
       double precision efczz(nefc)  !< [Input[ the COSMO zeta value for
@@ -490,6 +490,7 @@ c
       double precision ecos !< [Output] the energy contribution due to
                             !< the COSMO charges
       logical  status
+      logical  do_cosmo_file
 c
       logical  dbug,more,out,direct,noall,all,elec,nucl,iefc_done
       character*16 at_tag
@@ -1219,52 +1220,55 @@ c
 c
 c     ----- printing cosmo file for COSMO-RS/SAC codes -----
 c
-      cosmo_file_in = "cosmo"
-      call util_file_name(cosmo_file_in,.false.,.false.,cosmo_file) 
-      call util_file_name_resolve(cosmo_file,.false.)
-      if(ga_nodeid().eq.0) then
-        if(.not.util_io_unit(80,90,fn))
-     +     call errquit('cannot get free unit', 0,
-     +       0)
+      if (.not. rtdb_get(rtdb,'cosmo:do_cosmo_file',mt_log,1,
+     $     do_cosmo_file)) do_cosmo_file=.false.
+      if (do_cosmo_file) then
+         cosmo_file_in = "cosmo"
+         call util_file_name(cosmo_file_in,.false.,.false.,cosmo_file)
+         call util_file_name_resolve(cosmo_file,.false.)
+         if(ga_nodeid().eq.0) then
+         if(.not.util_io_unit(80,90,fn))
+     +     call errquit('cannot get free unit', 0, 0)
 c
-        open(unit=fn,form="formatted",file=cosmo_file)
-        if (dbug) then
-          write(*,*) "printing cosmo file for COSMO-RS/SAC",
+         open(unit=fn,form="formatted",file=cosmo_file)
+         if (dbug) then
+            write(*,*) "printing cosmo file for COSMO-RS/SAC",
      +     cosmo_file
-        end if
-        write(fn,"(a)") "$cosmo"
-        write(fn,*) " epsilon= ", dielec
-        write(fn,*) " rsolv= ", rsolv
+         end if
+         write(fn,"(a)") "$cosmo"
+         write(fn,*) " epsilon= ", dielec
+         write(fn,*) " rsolv= ", rsolv
 c
-        write(fn,"(a)") "$cosmo_data"
-        write(fn,*) " nps= ", nefc
-        write(fn,*) " area= ", srfmol/(cau2ang**2)
-        write(fn,*) " volume= ", volmol/(cau2ang**3)
+         write(fn,"(a)") "$cosmo_data"
+         write(fn,*) " nps= ", nefc
+         write(fn,*) " area= ", srfmol/(cau2ang**2)
+         write(fn,*) " volume= ", volmol/(cau2ang**3)
 c
-        write(fn,"(a)") "$coord_rad"
-        write(fn,"(a)") "#atom x y z [a.u.] element radius [A]"
-        do i = 1,nat
-          write(fn,*) i,
+         write(fn,"(a)") "$coord_rad"
+         write(fn,"(a)") "#atom x y z [a.u.] element radius [A]"
+         do i = 1,nat
+            write(fn,*) i,
      +     dbl_mb(k_xyzpt  +3*(i-1)),
      +     dbl_mb(k_xyzpt+1+3*(i-1)),
      +     dbl_mb(k_xyzpt+2+3*(i-1)),
      +     int(dbl_mb(k_zanpt+i-1)),
      +     ratm(i)*cau2ang
-        enddo ! i
+         enddo ! i
 c        
-        write(fn,"(a)") "$segment_information"
-        write(fn,"(a)") "# n atom x y z charge area charge/area"
-        write(fn,"(a)") "# coordinates [a.u.], area [A2]"
-        do ief=1,nefc
-           efcsi = efcs(ief)*(cau2ang**2)
-           if (do_cosmo_model.eq.DO_COSMO_YK)
-     +      efcsi = msrf(ief)*(cau2ang**2)
+         write(fn,"(a)") "$segment_information"
+         write(fn,"(a)") "# n atom x y z charge area charge/area"
+         write(fn,"(a)") "# coordinates [a.u.], area [A2]"
+         do ief=1,nefc
+            efcsi = efcs(ief)*(cau2ang**2)
+            if (do_cosmo_model.eq.DO_COSMO_YK)
+     +        efcsi = msrf(ief)*(cau2ang**2)
 
-           write(fn,*) ief, efciat(ief),
+            write(fn,*) ief, efciat(ief),
      +      efcc(1,ief), efcc(2,ief), efcc(3,ief),
      +      efcz(ief), efcsi, efcz(ief)/efcsi
-        end do
-        close(fn)
+         end do
+         close(fn)
+         end if
       end if
 c
 c MN solvation models -->