Index: INCL.H
===================================================================
--- INCL.H	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ INCL.H	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -2,4 +2,5 @@
       
       character*255 version
+      character*255 logString
 
       integer MAX_REPLICA, MAX_PROC, mxml, mxrs, mxat, mxbd, mxath
Index: INCP.H
===================================================================
--- INCP.H	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ INCP.H	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -44,12 +44,12 @@
 !      ir=0
 !      do i=1,nchp
-!        write(*,*) ' ===== chain |',chnp(i),'|'
+!        write (logString, *) ' ===== chain |',chnp(i),'|'
 !        do j=1,nchrsp(i)
 !          ir=ir+1
-!          write(*,*) ' ----- ',rsidp(ir),' ',rsnmp(ir),' ',nrsatp(ir)
+!          write (logString, *) ' ----- ',rsidp(ir),' ',rsnmp(ir),' ',nrsatp(ir)
 !          k1=irsatp(ir)
 !          k2=k1+nrsatp(ir)-1
 !          do k=k1,k2
-!            write(*,*) ' ',noatp(k),' ',atnmp(k),' ',(xyzp(l,k),l=1,3)
+!            write (logString, *) ' ',noatp(k),' ',atnmp(k),' ',(xyzp(l,k),l=1,3)
 !          enddo
 !        enddo
Index: addend.f
===================================================================
--- addend.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ addend.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -42,6 +42,7 @@
      &.or.grn(:3).eq.'nme'.or.grc(:3).eq.'nme') then
 
-        write(*,'(2a)') ' addend> N-Acetyl (ace) or N-Methylamide (nme)'
-     &   ,' should be put in SEQUENCE file, not added as end groups'
+        write (logString, '(2a)') 
+     &   ' addend> N-Acetyl (ace) or N-Methylamide (nme)',
+     &   ' should be put in SEQUENCE file, not added as end groups'
 
         stop
@@ -59,5 +60,5 @@
               sbrs='nh2+'
             else
-              write (*,'(2a)') ' addend> ',
+              write (logString, '(2a)') ' addend> ',
      &         ' No N-terminal Hyp possible with ECEPP/3 dataset'
               stop
@@ -87,5 +88,5 @@
         else
 
-          write(*,'(2a)') ' addend> Can add only ',
+          write (logString, '(2a)') ' addend> Can add only ',
      &     'nh2 or nh3+ to N-terminus'
           stop
@@ -98,5 +99,5 @@
       else ! ace
 
-        write(*,'(2a)') ' addend> Acetyl group',
+        write (logString, '(2a)') ' addend> Acetyl group',
      &     ' at N-terminus not modified'
       endif
@@ -128,5 +129,5 @@
         else
 
-          write(*,'(2a)') ' addend> Can add only ',
+          write (logString, '(2a)') ' addend> Can add only ',
      &     'cooh or coo- to C-terminus'
           stop
@@ -139,5 +140,5 @@
       else  ! N'-methylamide
 
-        write(*,'(2a)') ' addend> N-Methylamide',
+        write (logString, '(2a)') ' addend> N-Methylamide',
      &     ' at C-terminus not modified'
 
@@ -149,7 +150,9 @@
         cg = cg + cgat(i)
       enddo
-      if (abs(cg).gt.1.d-5) write(*,'(a,i2,a,f7.3,/)')
-     &        ' addend> Net charge of molecule #'
+      if (abs(cg).gt.1.d-5) then
+         write (logString, '(a,i2,a,f7.3)') 
+     &      ' addend> Net charge of molecule #'
      &        ,nml,': ',cg
+      endif
 
       return
@@ -204,6 +207,7 @@
         endif
       enddo
-      write (*,'(4a,i4,a,i4)') ' rplgrp> cannot find atom >',rpat,
-     &'< to be replaced in residue ',seq(nrs),nrs,' of molecule ',nml
+      write (logString, '(4a,i4,a,i4)') ' rplgrp> cannot find atom >',
+     &   rpat,
+     &   '< to be replaced in residue ',seq(nrs),nrs,' of molecule ',nml
       stop
 
@@ -234,5 +238,5 @@
             if (iowat(ibd(i)).ne.nfirp) then
               if (ibdrg.ne.0) then
-                write (*,'(2a,i3)') 
+                write (logString, '(2a,i3)') 
      &             ' rplgrp> Can handle only simple ring at 1st',
      &             ' atom of molecule #',nml
@@ -261,5 +265,5 @@
           goto 11
         else
-          write (*,'(4a,i4,a,i4)') 
+          write (logString, '(4a,i4,a,i4)') 
      &      ' rplgrp> Cannot replace BACKBONE atom ',rpat,
      &      ' of residue ',seq(nrs),nrs,' in molecule #',nml
@@ -275,5 +279,5 @@
 ! _______________________________ get data for substituent atoms
     3 if (iopfil(lunlib,reslib,'old','formatted').le.izero) then
-        write (*,'(a,/,a,i3,2a)') 
+        write (logString, '(a,/,a,i3,2a)') 
      &    ' rplgrp> ERROR opening library of residues:',
      &    ' LUN=',lunlib,' FILE=',reslib(1:iendst(reslib))
@@ -289,5 +293,5 @@
         endif
       enddo
-      write (*,'(4a)') ' rplgrp> Cannot find atom >',rpat,
+      write (logString, '(4a)') ' rplgrp> Cannot find atom >',rpat,
      &'< in substituent residue ',sbrs
       stop
@@ -372,5 +376,6 @@
         else  ! more atoms
           if ((ilaat+nsh).gt.mxat) then
-            write (*,'(a,i5)') ' rplgrp> number of atoms > ',mxat
+            write (logString, '(a,i5)') ' rplgrp> number of atoms > ',
+     &         mxat
             stop
           endif
@@ -466,5 +471,5 @@
           nb=nb+1
           if (nb.gt.mxbd) then
-            write (*,'(6a,/,2a,3(i4,a))')
+            write (logString, '(6a,/,2a,3(i4,a))')
      &      ' rplgrp> Cannot add atoms following ',rpat,
      &      ' from group ',sbrs,' to atom ',rpat,
@@ -593,5 +598,6 @@
         else
           if ((ilavr+jsh).gt.mxvr) then
-            write (*,'(a,i5)') ' rplgrp> number of variables > ',mxvr
+            write (logString, '(a,i5)') 
+     &         ' rplgrp> number of variables > ',mxvr
             stop
           endif
@@ -628,10 +634,10 @@
       return
 ! __________________________________________ Errors
-   10 write (*,'(3a,/,2a,i4,a,i4,/,2a)') 
+   10 write (logString, '(3a,/,2a,i4,a,i4,/,2a)') 
      &   ' rplgrp> Cannot replace atom(s) following ',rpat,
      &   ' from INSIDE a ring','    in residue: ',seq(nrs),nrs,
      &   ' in molecule #',nml,' or in substitute: ',sbrs
       stop
-   11 write (*,'(4a,i4,a,i4,/,a)') 
+   11 write (logString, '(4a,i4,a,i4,/,a)') 
      &   ' rplgrp> Cannot replace atom(s) following ',rpat,
      &   ' of residue ',seq(nrs),nrs,' in molecule #',nml,
@@ -686,5 +692,6 @@
           endif
         else
-          write (*,*) ' redchg> dont know which end goup is present'
+          write(logString,*)
+     &                  ' redchg> dont know which end goup is present'
           stop
         endif
@@ -698,5 +705,5 @@
 
         if (iopfil(lunchg,chgfil,'old','formatted').le.izero) then
-          write (*,'(a,/,a,i3,2a)') 
+          write (logString, '(a,/,a,i3,2a)') 
      &      ' redchg> ERROR opening library of charges:',
      &      ' LUN=',lunchg,' FILE=',chgfil(1:iendst(chgfil))
@@ -722,5 +729,6 @@
                 endif
               enddo
-              write (*,'(6a)') ' redchg> Cannot find atom: ',atnm,
+              write (logString, '(6a)') ' redchg> Cannot find atom: ',
+     &                        atnm,
      &                        ' for entry: ',cgty,' in library: ',
      &                        chgfil(1:iendst(chgfil))
@@ -729,5 +737,5 @@
             return
           else
-            write (*,'(4a)')
+            write (logString, '(4a)')
      &       ' redchg> must increase MXATH to read data for entry: ',
      &        cgty,' in library: ',chgfil(1:iendst(chgfil))
@@ -737,5 +745,5 @@
         endif
         goto 1
-    3   write (*,'(4a)')
+    3   write (logString, '(4a)')
      &   ' redchg> Cannot find entry: ',cgty,' in library: ',
      &      chgfil(1:iendst(chgfil))
@@ -746,5 +754,5 @@
 
         if (iopfil(lunlib,reslib,'old','formatted').le.izero) then
-          write (*,'(a,/,a,i3,2a)') 
+          write (logString, '(a,/,a,i3,2a)') 
      &      ' redchg> ERROR opening library of residues:',
      &      ' LUN=',lunlib,' FILE=',reslib(1:iendst(reslib))
@@ -769,5 +777,6 @@
                 endif
               enddo
-              write (*,'(6a)') ' redchg> Cannot find atom: ',atnm,
+              write (logString, '(6a)') ' redchg> Cannot find atom: ',
+     &                        atnm,
      &                        ' for entry: ',cgty,' in library: ',
      &                        reslib(1:iendst(reslib))
@@ -776,5 +785,5 @@
             return
           else
-            write (*,'(4a)')
+            write (logString, '(4a)')
      &       ' redchg> must increase MXATH to read data for entry: ',
      &        cgty,' in library: ',reslib(1:iendst(reslib))
@@ -784,5 +793,5 @@
         endif
         goto 4
-    6   write (*,'(4a)')
+    6   write (logString, '(4a)')
      &   ' redchg> Cannot find entry: ',cgty,' in library: ',
      &      reslib(1:iendst(reslib))
@@ -792,5 +801,5 @@
       endif
 
-   10 write (*,'(4a)')
+   10 write (logString, '(4a)')
      &    ' redchg> Do not have charges for N/C-terminal residue ',
      &    res,' modified with group :',sbrs
Index: anneal.f
===================================================================
--- anneal.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ anneal.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -84,5 +84,6 @@
 
       eol=energy()
-      write (*,'(a,e12.5,/)')  'energy of start configuration: ',eol
+      write (logString, '(a,e12.5,/)')  
+     &   'energy of start configuration: ',eol
 
 ! Write start configuration in pdb-format into file
@@ -94,5 +95,6 @@
          call metropolis(eol,acz,can_weight)
       end do
-      write(*,*) 'Energy after  equilibration:',eol
+      write (logString, *) 
+     &   'Energy after  equilibration:', eol
 
 !======================Simulation by simulated annealing
@@ -131,15 +133,15 @@
 
       acz = acz/dble(nsw*nvr)
-      write(*,*) 'acceptance rate:',acz
-      write(*,*)
+      write (logString, *) 'acceptance rate:',acz
+      write (logString, *)
 ! ------------ Output Dihedreals of final configuration
-      write(*,*) 'last energy',eol
+      write (logString, *) 'last energy',eol
       call outvar(0,' ')
 !     Output final conformation as pdb-file
       call outpdb(0,"final.pdb")
-      write(*,*)
+      write (logString, *)
 
 ! ------------ Output Dihedreals of conformation with lowest energy
-      write(*,*) 'lowest energy ever found:',nemin,ymin
+      write (logString, *) 'lowest energy ever found:',nemin,ymin
       close(14)
 ! =====================
Index: bgs.f
===================================================================
--- bgs.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ bgs.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -327,7 +327,7 @@
 !          enw=energy()
 !          if (abs(enw-eol1).gt.0.000001) then
-!             write(*,*) 'rejected bgs move: energy change :',eol1,enw
+!             write (logString, *) 'rejected bgs move: energy change :',eol1,enw
 !          endif
-!         write(*,*) 'rejected bgs move: ',eol1,enw,wfw,wbw,rd
+!         write (logString, *) 'rejected bgs move: ',eol1,enw,wfw,wbw,rd
          bgs=0
       endif
Index: bldmol.f
===================================================================
--- bldmol.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ bldmol.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -134,5 +134,5 @@
 
           if (ix.ne.0) then
-            write (*,'(2a,i3)') 
+            write (logString, '(2a,i3)') 
      &         ' fnd3ba> Can handle only simple ring at 1st',
      &         ' atom of molecule #',nml
@@ -166,5 +166,5 @@
       endif
 
-      write (*,'(4a,i4,a,i4)') 
+      write (logString, '(4a,i4,a,i4)') 
      &   ' fnd3ba> Cannot find backbone atom following ',nmat(i2),
      &   ' of residue ',seq(irs),irs,' in molecule #',nml
Index: canon.f
===================================================================
--- canon.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ canon.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -69,5 +69,5 @@
       eol = energy()
       e_min = eol
-      write (*,'(a,e12.5,/)')  'energy of start configuration:',eol
+      write (logString, '(a,e12.5,/)')  'energy of start configuration:',eol
 
 ! Write start configuration in pdb-format into file
@@ -79,5 +79,5 @@
          call metropolis(eol,acz,can_weight)
       end do
-      write(*,*) 'Energy after equilibration:',eol
+      write (logString, *) 'Energy after equilibration:',eol
 
 !======================Simulation in canonical ensemble
@@ -86,6 +86,6 @@
         call metropolis(eol,acz,can_weight)
         if (eol.lt.e_min) then
-           write (*,*) "New minimum energy:", eol, "t=", nsw
-           write (*,*) eyel,eyhb,eyvr,eysl
+           write (logString, *) "New minimum energy:", eol, "t=", nsw
+           write (logString, *) eyel,eyhb,eyvr,eysl
            e_min = eol
            call outpdb(0, "minen.pdb")
@@ -116,8 +116,8 @@
 
       acz = acz/dble(nsw*nvr)
-      write(*,*) 'acceptance rate:',acz
-      write(*,*)
+      write (logString, *) 'acceptance rate:',acz
+      write (logString, *)
 ! ------------ Output Dihedreals of final configuration
-      write(*,*) 'last energy',eol
+      write (logString, *) 'last energy',eol
       call outvar(0,'lastconf.var')
 !     Output final conformation as pdb-file
Index: cnteny.f
===================================================================
--- cnteny.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ cnteny.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -37,5 +37,5 @@
 
       if (ntlvr.eq.0) then
-        write (*,'(a,i4)')
+        write (logString, '(a,i4)')
      &           ' cnteny> No variables defined in molecule #',nml
         return
@@ -184,6 +184,6 @@
           nbc=nbc+1
           ir=nursat(i)
-          write(*,'(1x,i4,1x,a4,1x,a4,a2,e11.4)') ir,seq(ir),nmat(i),
-     &                                            ': ',ey
+          write (logString, '(1x,i4,1x,a4,1x,a4,a2,e11.4)') ir,seq(ir),
+     &       nmat(i), ': ',ey
         endif
       enddo
Index: contacts.f
===================================================================
--- contacts.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ contacts.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -161,5 +161,5 @@
            do nr_j=nr_i+3,irsml2(nmol) ! Over res. j 
 
-!             write(*,'(2i3)'),nr_i,nr_j
+!             write (logString, '(2i3)'),nr_i,nr_j
 
               ic=0
@@ -173,5 +173,5 @@
               if(sqrt(rij2).lt.rcut) ic=1
 
-!             write(*,'(2i3)'),nr_i,nr_j
+!             write (logString, '(2i3)'),nr_i,nr_j
 
               ijcont(nr_i,nr_j)=ic
Index: dihedr.f
===================================================================
--- dihedr.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ dihedr.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -57,11 +57,12 @@
         return
       else
-        write (*,'(a,4i5)')' dihedr> Error in coordinates of atoms #: '
+        write (logString, '(a,4i5)')
+     &     ' dihedr> Error in coordinates of atoms #: '
      &                     ,i1,i2,i3,i4
 
-        write (*,*) 'stored coordinates are xvals :',
+        write (logString, *) 'stored coordinates are xvals :',
      &       xat(i1),xat(i2),xat(i3),xat(i4) 
-        write (*,*) 'yvals:', yat(i1),yat(i2),yat(i3),yat(i4) 
-        write (*,*) 'zvals:', zat(i1),zat(i2),zat(i3),zat(i4) 
+        write (logString, *) 'yvals:', yat(i1),yat(i2),yat(i3),yat(i4) 
+        write (logString, *) 'zvals:', zat(i1),zat(i2),zat(i3),zat(i4) 
         call outvar(0,'crash.var')        
         stop
@@ -110,10 +111,11 @@
 
       else
-        write (*,'(a,3i5)')' valang> Error in coordinates of atoms #: '
+        write (logString, '(a,3i5)')
+     &   ' valang> Error in coordinates of atoms #: '
      &                     ,i1,i2,i3
-        write (*,*) 'stored coordinates are xvals :',
+        write (logString, *) 'stored coordinates are xvals :',
      &       xat(i1),xat(i2),xat(i3) 
-        write (*,*) 'yvals:', yat(i1),yat(i2),yat(i3) 
-        write (*,*) 'zvals:', zat(i1),zat(i2),zat(i3) 
+        write (logString, *) 'yvals:', yat(i1),yat(i2),yat(i3) 
+        write (logString, *) 'zvals:', zat(i1),zat(i2),zat(i3) 
         call outvar(0,'crash.var') 
         stop
Index: enyflx.f
===================================================================
--- enyflx.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ enyflx.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -33,5 +33,5 @@
       ntlvr=nvrml(nml)
       if (ntlvr.eq.0) then
-        write (*,'(a,i4)')
+        write (logString, '(a,i4)')
      &           ' enyflx> No variables defined in molecule #',nml
         return
Index: enyshe.f
===================================================================
--- enyshe.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ enyshe.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -41,5 +41,5 @@
 
       if (ntlvr.eq.0) then
-        write (*,'(a,i4)')
+        write (logString, '(a,i4)')
      &           ' enyshe> No variables defined in molecule #',nml
         return
Index: enyshe_p.f
===================================================================
--- enyshe_p.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ enyshe_p.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -46,5 +46,5 @@
 
       if (ntlvr.eq.0) then
-         write (*,'(a,i4)')
+         write (logString, '(a,i4)')
      &        ' enyshe> No variables defined in molecule #',nml
          return
Index: enysol.f
===================================================================
--- enysol.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ enysol.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -30,4 +30,5 @@
 !     functions
       integer nursat
+
 
       integer numbox, inbox, indsort, look, i, ii, ia, ib, ibox, icount
@@ -157,10 +158,11 @@
         numbox(j)=nboxj
         if (nboxj.gt.mxbox) then
-         write(*,'(a)') 'enysol> bad mxboxe-2'
-         write(*,*) 'diagnostics ...'
-         write(*,*) 'atom ',j
-         write(*,*) 'position ',xat(j),yat(j),zat(j)
-         write(*,*) 'box indices ',mx,my,mz
-         write(*,*) 'resulting boxindex and limit ',nboxj,mxbox
+         write (logString, '(a)') 'enysol> bad mxboxe-2'
+         write (logString, *) 'diagnostics ...'
+         write (logString, *) 'atom ',j
+         write (logString, *) 'position ',xat(j),yat(j),zat(j)
+         write (logString, *) 'box indices ',mx,my,mz
+         write (logString, *) 'resulting boxindex and limit ',
+     &      nboxj,mxbox
 
          stop
@@ -350,5 +352,5 @@
 
       read(lin(1:5),'(i5)') npnt
-!       write(*,'(a,i5)') 'the number of points---->',npnt
+!       write (logString, '(a,i5)') 'the number of points---->',npnt
 
 !    Read the surface points
Index: enysol_p.f
===================================================================
--- enysol_p.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ enysol_p.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -159,10 +159,11 @@
         numbox(j)=nboxj
         if (nboxj.gt.mxbox) then
-         write(*,'(a)') 'enysol> bad mxboxe-2'
-         write(*,*) 'diagnostics ...'
-         write(*,*) 'atom ',j
-         write(*,*) 'position ',xat(j),yat(j),zat(j)
-         write(*,*) 'box indices ',mx,my,mz
-         write(*,*) 'resulting boxindex and limit ',nboxj,mxbox
+         write (logString, '(a)') 'enysol> bad mxboxe-2'
+         write (logString, *) 'diagnostics ...'
+         write (logString, *) 'atom ',j
+         write (logString, *) 'position ',xat(j),yat(j),zat(j)
+         write (logString, *) 'box indices ',mx,my,mz
+         write (logString, *) 'resulting boxindex and limit ',nboxj,
+     &      mxbox
 
          stop
@@ -195,6 +196,6 @@
       if (myrank.eq.(no - 1)) iboxmax = ncbox
       if (myrank.eq.-1) then
-         write (*,*) 'enysol> Summary:', enysolct, ncbox, boxpp,
-     &               ndx, ndy, ndz
+         write (logString, *) 'enysol> Summary:', enysolct, ncbox, 
+     &               boxpp, ndx, ndy, ndz
       endif
       do ibox = iboxmin + 1, iboxmax + 1
@@ -333,5 +334,5 @@
       call MPI_ALLREDUCE(eysl, eyslsum, 1, MPI_DOUBLE_PRECISION,
      &      MPI_SUM,my_mpi_comm, ierror)
-!       write(*,*) 'enysol>', myrank, eysl, eyslsum
+!       write (logString, *) 'enysol>', myrank, eysl, eyslsum
       tsurfres = surfres
       call MPI_ALLREDUCE(tsurfres, surfres, mxrs, MPI_DOUBLE_PRECISION,
@@ -341,5 +342,5 @@
       endwtime = MPI_Wtime()
       if (myrank.le.-1) then
-         write (*,*) 'enysol>',myrank,enysolct,
+         write (logString, *) 'enysol>',myrank,enysolct,
      &               iboxmin + 1, iboxmax + 1,
      &               endwtime - startwtime, "s"
@@ -376,5 +377,5 @@
 
       read(lin(1:5),'(i5)') npnt
-!        write(*,'(a,i5)') 'the number of points---->',npnt
+!        write (logString, '(a,i5)') 'the number of points---->',npnt
 
 !    Read the surface points
Index: esolan.f
===================================================================
--- esolan.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ esolan.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -130,6 +130,7 @@
      &   +(zold(j)-zold(i))**2! Distance between atom i and j
        if(ddp.lt.1e-10) then
-         write(*,*)'ERROR in data: centres of two atoms coincide!'
-         write(*,*)i,j,xold(i),yold(i),zold(i),rvdw(i),
+         write (logString, *)
+     &      'ERROR in data: centres of two atoms coincide!'
+         write (logString, *)i,j,xold(i),yold(i),zold(i),rvdw(i),
      &                   xold(j),yold(j),zold(j),rvdw(j)
          stop 'Centres of atoms coincide!'
Index: getmol.f
===================================================================
--- getmol.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ getmol.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -49,5 +49,5 @@
 
       if (iopfil(lunlib,reslib,'old','formatted').le.izero) then
-        write (*,'(a,/,a,i3,2a)') 
+        write (logString, '(a,/,a,i3,2a)') 
      &    ' getmol> ERROR opening library of residues:',
      &    ' LUN=',lunlib,' FILE=',reslib(1:iendst(reslib))
@@ -74,10 +74,10 @@
 
         if (res(:3).eq.'nme'.and.nrs.ne.ilars) then
-          write (*,'(3a)') ' getmol> residue >',res,
+          write (logString, '(3a)') ' getmol> residue >',res,
      &                     '< allowed at C-terminus only !'
           close(lunlib)
           stop
         elseif (res(:3).eq.'ace'.and.nrs.ne.ifirs) then
-          write (*,'(3a)') ' getmol> residue >',res,
+          write (logString, '(3a)') ' getmol> residue >',res,
      &                     '< allowed at N-terminus only !'
           close(lunlib)
@@ -88,10 +88,11 @@
 
         if ((nat+ntlat).gt.mxat) then
-          write (*,'(a,i5)') ' getmol> number of atoms > ',mxat
+          write (logString, '(a,i5)') ' getmol> number of atoms > ',mxat
           close(lunlib)
           stop
         endif
         if ((nvr+ntlvr).gt.mxvr) then
-          write (*,'(a,i5)') ' getmol> number of variables > ',mxvr
+          write (logString, '(a,i5)') ' getmol> number of variables > '
+     &       ,mxvr
           close(lunlib)
           stop
@@ -139,5 +140,5 @@
               nbd=nbdat(nh)
               if (nbd.eq.mxbd) then
-                write(*,'(a,i2,a,i4,2a,i4,a)') 
+                write (logString, '(a,i2,a,i4,2a,i4,a)') 
      &           ' getmol> need ',(mxbd+2),
      &           'th bond to connect residues ',
@@ -154,5 +155,5 @@
 
                   if (t.eq.0.0) then
-                    write (*,'(3a,/,2a)') 
+                    write (logString, '(3a,/,2a)') 
      &               ' getmol> DIHEDRAL for atom ',nmat(nj),
      &               ' should be PHASE angle with respect to atom ',
@@ -312,5 +313,5 @@
 
         if (nat.gt.mxath) then
-          write (*,'(a,i5)') ' redres> number of atoms > ',mxath
+          write (logString, '(a,i5)') ' redres> number of atoms > ',mxath
           close(lunlib)
           stop
@@ -340,8 +341,9 @@
           if (iexcp.eq.0.and.i.le.jow) then
             if (i.eq.jow) then
-              write (*,'(5a)') ' redres> atom ',nmath(i),' of ',
+              write(logString,'(5a)') ' redres> atom ',nmath(i),' of ',
      &                          resl,' cannot preceed itself '
             else
-              write (*,'(5a,i4)') ' redres> atom ',nmath(i),' of ',
+              write (logString, '(5a,i4)') 
+     &                    ' redres> atom ',nmath(i),' of ',
      &                    resl,' should be placed AFTER atom #',jow
             endif
@@ -394,5 +396,6 @@
 
               if (nvrr.gt.mxvrh) then
-                write (*,'(a,i5)') ' redres> number of variables > ',
+                write (logString, '(a,i5)') 
+     &           ' redres> number of variables > ',
      &                             mxvrh
                 close(lunlib)
@@ -423,28 +426,31 @@
 
 ! ____________________________________________________________ ERRORS
-    2 write (*,'(4a)') ' redres> residue >',resl,'< NOT FOUND in ',
+    2 write (logString, '(4a)') ' redres> residue >',resl,
+     &      '< NOT FOUND in ',
      &reslib(1:iendst(reslib))
       close(lunlib)
       stop
 
-    3 write (*,'(a,i4,2a)') ' redres> ERROR reading line No. ',nln,
-     &' in ',reslib(1:iendst(reslib))
+    3 write (logString, '(a,i4,2a)') ' redres> ERROR reading line No. ',
+     &   nln,
+     &   ' in ',reslib(1:iendst(reslib))
       close(lunlib)
       stop
 
-    4 write (*,'(4a)') ' redres> Incorrect order of bonds for atom ',
+    4 write (logString, '(4a)') 
+     &   ' redres> Incorrect order of bonds for atom ',
      &                      nmath(i),' of ',resl
 
-    5 write (*,'(8x,2a)') '... must correct ',
+    5 write (logString, '(8x,2a)') '... must correct ',
      &                      reslib(1:iendst(reslib))
       close(lunlib)
       stop
 
-    6 write (*,'(a,i2,4a)') ' redres> unknown type :',ity,
+    6 write (logString, '(a,i2,4a)') ' redres> unknown type :',ity,
      &                   ': for atom ',nmath(i),' in residue ',resl
       close(lunlib)
       stop
 
-    7 write (*,'(a,i2,4a)') ' redres> unknown class :',ic,
+    7 write (logString, '(a,i2,4a)') ' redres> unknown class :',ic,
      &                   ': for variable ',nm(j),' in residue ',resl
       close(lunlib)
Index: gradient.f
===================================================================
--- gradient.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ gradient.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -57,5 +57,5 @@
         else
 
-          write(*,*)  'gradient> Set itysol < 0'
+          write (logString, *)  'gradient> Set itysol < 0'
           stop
         endif
Index: hbond.f
===================================================================
--- hbond.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ hbond.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -59,5 +59,5 @@
         ntlvr=nvrml(nml)
         if (ntlvr.eq.0) then
-          write (*,'(a,i4)')
+          write (logString, '(a,i4)')
      &           ' hbond> No variables defined in molecule #',nml
           return
@@ -108,5 +108,6 @@
                     ixtyhb(ntyhb)=ixhb
                   else
-                    write(*,*) ' hbond> increase parameter MXTYHB'
+                    write (logString, *) 
+     &                    ' hbond> increase parameter MXTYHB'
                     stop
                   endif
@@ -139,5 +140,6 @@
                   ixtyhb(ntyhb)=ixhb
                 else
-                  write(*,*) ' hbond> increase parameter MXTYHB'
+                  write (logString, *) 
+     &               ' hbond> increase parameter MXTYHB'
                   stop
                 endif
@@ -158,5 +160,6 @@
 !     enddo
 
-      if (ipr.gt.0) write(*,'(1x,a,/)') ' hbond>  Hydrogen Bonds:'
+      if (ipr.gt.0) write (logString, '(1x,a,/)') 
+     &   ' hbond>  Hydrogen Bonds:'
 
       if (ntyhb.gt.0) then
@@ -287,8 +290,8 @@
       adab=valang(id,ia,ib)*crd
 
-      write(*,*) '  '
-      write(*,*) ' Dah: ',dah,' Dad: ',dad
-      write(*,*) ' Adha: ',adha,' Ahab: ',adha,' Adab: ',adab
-      write(*,*) '  '
+      write (logString, *) '  '
+      write (logString, *) ' Dah: ',dah,' Dad: ',dad
+      write (logString, *) ' Adha: ',adha,' Ahab: ',adha,' Adab: ',adab
+      write (logString, *) '  '
 
       return
Index: init_energy.f
===================================================================
--- init_energy.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ init_energy.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -66,6 +66,6 @@
 
 !----Initialize solvation part if necessary
-      write (*,*) 'init_energy: itysol = ',itysol
-      write(*,*) 'init_energy: esol_scaling = ',isolscl
+      write (logString, *) 'init_energy: itysol = ',itysol
+      write (logString, *) 'init_energy: esol_scaling = ',isolscl
 
       its = iabs(itysol)
@@ -82,5 +82,6 @@
 
         if (itysol.ne.0) then
-          write(*,'(a)') ' init_energy>  undefined solvent type !'
+          write (logString, '(a)') 
+     &      ' init_energy>  undefined solvent type !'
           stop
         endif
@@ -96,5 +97,6 @@
       return
 
- 10   write (*, '(a)') 'Cannot open the file with surface points'
+ 10   write (logString,  '(a)') 
+     &   'Cannot open the file with surface points'
       stop
 
Index: main.f
===================================================================
--- main.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ main.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -81,5 +81,5 @@
 
       ntlml = 0
-      write (*,*) 'Solvent: ', itysol
+      write (logString, *) 'Solvent: ', itysol
 !     Initialize random number generator.
       call sgrnd(31433)
@@ -118,8 +118,8 @@
        lrand = .true.
        E = energy()
-       write(*,*) E, eyel,eyvw,eyhb,eyvr
+       write (logString, *) E, eyel,eyvw,eyhb,eyvr
        call outpdb(1, "polyA.pdb")
 !      Canonical Monte Carlo
-!        call canon(nequi, nsweep, nmes, temp, lrand)
+       call canon(nequi, nsweep, nmes, temp, lrand)
 
 !      For simulated annealing uncomment the lines below
Index: main_bgl_p.f
===================================================================
--- main_bgl_p.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ main_bgl_p.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -57,5 +57,5 @@
       call mpi_init(ierr)
 !       call pmi_cart_comm_create(comm_cart,ierr)
-      write(*,*) "Initialized MPI. Now setting up communicators." 
+      write (logString, *) "Initialized MPI. Now setting up communicators." 
       call flush(6)
       ndims = 4
@@ -81,5 +81,5 @@
       call MPI_Cart_GET(comm_cart, ndims, dims, periods, coords, ierr)
 
-      write(*,*) ndims, dims, periods, coords
+      write (logString, *) ndims, dims, periods, coords
       call flush(6)
 !       call VTSetup()
@@ -192,5 +192,5 @@
      &      + (coords(3)/ldims(3))*nblock(1)*nblock(2)
 
-      write(*,*) myrank, color, ldims, nblock
+      write (logString, *) myrank, color, ldims, nblock
 
       call mpi_comm_split(comm_cart,color,myrank,local_comm,ierr)
@@ -216,5 +216,5 @@
 !          call mpi_group_range_incl(group_world, 1, proc_range, group(i)
 !      &                              ,error)
-!          write (*,*) "Assigning rank ", j, proc_range, 
+!          write (logString, *) "Assigning rank ", j, proc_range, 
 !      &               "to group", group(i)
 !          call flush(6)
@@ -227,5 +227,5 @@
 !              my_mpi_comm = comm(i)
 !              rep_id = i - 1
-!              write (*,*) rep_id, "has comm", my_mpi_comm
+!              write (logString, *) rep_id, "has comm", my_mpi_comm
 !              call flush(6)
 !          endif
@@ -233,5 +233,5 @@
 ! 
 ! c     Setup the communicator used for parallel tempering
-!       write (*,*) "PTGroup=", ranks(:num_replica)
+!       write (logString, *) "PTGroup=", ranks(:num_replica)
 !       call flush(6)
 !       call mpi_group_incl(group_world, num_replica, ranks, group_partem,
@@ -241,5 +241,5 @@
 ! 
 !       if (partem_comm.ne.MPI_COMM_NULL) then
-!          write (*,*) partem_comm,myrank, "is master for ", rep_id, "."
+!          write (logString, *) partem_comm,myrank, "is master for ", rep_id, "."
 !       endif
 
@@ -247,9 +247,9 @@
       call mpi_comm_size(my_mpi_comm,no,ierr)
       rep_id = color
-      write (*,*) "My new rank is ", myrank, "of", no
+      write (logString, *) "My new rank is ", myrank, "of", no
       call flush(6)
       if (myrank.eq.0) then 
          color = 1
-         write (*,*) 'My rank and color:', myrank, color
+         write (logString, *) 'My rank and color:', myrank, color
          call flush(6)
       else 
@@ -258,5 +258,5 @@
       call mpi_comm_split(comm_cart,color,0,partem_comm,ierr)
 
-!       write(*,*) "Finalizing MPI."
+!       write (logString, *) "Finalizing MPI."
 !       call flush(6)
 !       CALL mpi_finalize(ierr)
@@ -295,5 +295,5 @@
 
 !     ========================================  start of parallel tempering run
-      write (*,*) "There are ", no,
+      write (logString, *) "There are ", no,
      &            " processors available for ",rep_id
       call flush(6)
@@ -309,5 +309,5 @@
 
       if(my_pt_rank.eq.0) then
-         write(*,*) "time for simulation using ", num_proc,
+         write (logString, *) "time for simulation using ", num_proc,
      &        " processors =",  endwtime - startwtime, " seconds"
          call flush(6)
Index: main_p.f
===================================================================
--- main_p.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ main_p.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -151,5 +151,5 @@
          call mpi_group_range_incl(group_world, 1, proc_range, group(i)
      &                              ,error)
-         write (*,*) "Assigning rank ", j, proc_range,
+         write (logString, *) "Assigning rank ", j, proc_range,
      &               "to group", group(i)
          call flush(6)
@@ -162,5 +162,5 @@
              my_mpi_comm = comm(i)
              rep_id = i - 1
-             write (*,*) rep_id, "has comm", my_mpi_comm
+             write (logString, *) rep_id, "has comm", my_mpi_comm
              call flush(6)
          endif
@@ -168,5 +168,5 @@
 
 !     Setup the communicator used for parallel tempering
-      write (*,*) "PTGroup=", ranks(:num_replica)
+      write (logString, *) "PTGroup=", ranks(:num_replica)
       call flush(6)
       call mpi_group_incl(group_world, num_replica, ranks, group_partem,
@@ -176,5 +176,6 @@
 
       if (partem_comm.ne.MPI_COMM_NULL) then
-         write (*,*) partem_comm,myrank, "is master for ", rep_id, "."
+         write (logString, *) partem_comm,myrank, "is master for ", 
+     &      rep_id, "."
       endif
 
@@ -182,5 +183,5 @@
       call mpi_comm_size(my_mpi_comm,no,ierr)
 
-      write (*,*) "My new rank is ", myrank, "of", no
+      write (logString, *) "My new rank is ", myrank, "of", no
       call flush(6)
 ! = Done setting up communicators =====================================
@@ -220,5 +221,5 @@
 
 !     ========================================  start of parallel tempering run
-      write (*,*) "There are ", no,
+      write (logString, *) "There are ", no,
      &            " processors available for ",rep_id
       call flush(6)
@@ -234,5 +235,5 @@
 
       if(my_pt_rank.eq.0) then
-         write(*,*) "time for simulation using ", num_proc,
+         write (logString, *) "time for simulation using ", num_proc,
      &        " processors =",  endwtime - startwtime, " seconds"
          call flush(6)
Index: metropolis.f
===================================================================
--- metropolis.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ metropolis.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -172,5 +172,6 @@
       enw = energy()
       if(abs(eol-enw).gt.0.000001)  then
-         write(*,*) 'metropolis: eol,enw,difference:', eol, enw, eol-enw
+         write (logString, *) 'metropolis: eol,enw,difference:', eol, 
+     &      enw, eol-enw
          if (eol.lt.100000) then 
             stop 'metropolis: eol and enw difference unacceptable'
Index: mincjg.f
===================================================================
--- mincjg.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ mincjg.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -57,5 +57,6 @@
      &          d(mxn),xa(mxn),ga(mxn),dt(mxn),yt(mxn),gt(mxn)
 
-
+      character(255) logString
+      
       ier = 0
       iter = 0
@@ -93,5 +94,5 @@
 
       if ( gdit .ge. 0.d0 ) then
-        write(*,*) ' mincjg>  search direction is uphill'
+        write (logString, *) ' mincjg>  search direction is uphill'
         ier = 3
         goto 6
@@ -162,5 +163,6 @@
 
           ier=2
-          write(*,*) ' mincjg>  too small step in search direction'
+          write (logString, *) 
+     &       ' mincjg>  too small step in search direction'
         endif
 
@@ -230,5 +232,6 @@
       if ( iter .ge. (iterfm + MXFCON) ) then
         ier = 4
-        write(*,*) ' mincjg>  line search failed to reduce function'
+        write (logString, *) 
+     &    ' mincjg>  line search failed to reduce function'
         return
       endif
@@ -304,5 +307,5 @@
       call move(nfun,n,f,x,g)
 
-      write(*,*) ' mincjg> ier = ',ier
+      write (logString, *) ' mincjg> ier = ',ier
 
       return
Index: minim.f
===================================================================
--- minim.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ minim.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -58,5 +58,5 @@
 
       if (msv.gt.msvmx) then
-        write (*,'(a,i5)') ' minreg> parameter MSV > ',msvmx
+        write (logString, '(a,i5)') ' minreg> parameter MSV > ',msvmx
         stop
       endif
@@ -66,10 +66,10 @@
       call gradient()
 
-      write(*,'(/,a,/)') ' Energy BEFORE minimization:' 
+      write (logString, '(/,a,/)') ' Energy BEFORE minimization:' 
 
       if (ireg.eq.0) then
 
-        write (*,'(a,e12.5,/,3(a,e11.4),/,2(a,e11.4),/)') ' Total: ',
-     &    eysm,
+        write (logString, '(a,e12.5,/,3(a,e11.4),/,2(a,e11.4),/)') 
+     &    ' Total: ', eysm,
      &    '   Coulomb: ',eyel,' Lennard-Jones: ',eyvw,' HB: ',eyhb,
      &    '   Variables: ',eyvr,'  Solvatation: ',eysl
@@ -78,6 +78,6 @@
 
 
-        write (*,'(a,e12.5,/,3(a,e11.4),/,3(a,e11.4),/)') ' Total: ',
-     &    wtey*eysm + wtrg*eyrg,
+        write (logString, '(a,e12.5,/,3(a,e11.4),/,3(a,e11.4),/)') 
+     &    ' Total: ', wtey*eysm + wtrg*eyrg,
      &    '   Coulomb: ',eyel,' Lennard-Jones: ',eyvw,' HB: ',eyhb,
      &    '   Variables: ',eyvr,'  Solvatation: ',eysl,
@@ -148,6 +148,6 @@
 !          gdrg2 = max(acc,gdrg2)
 !          wtrg = wtrg * sqrt(gdey2/gdrg2)
-!          write(*,*)  ' -->    Wt_energy = ',wtey,'  Wt_regul. = ',wtrg
-!          write(*,*)  '  '
+!          write (logString, *)  ' -->    Wt_energy = ',wtey,'  Wt_regul. = ',wtrg
+!          write (logString, *)  '  '
 !        endif
 
@@ -184,11 +184,11 @@
 
       if (nop.lt.mxop) then
-        write (*,'(a)') ' ---- CONVERGENCE ----'
+        write (logString, '(a)') ' ---- CONVERGENCE ----'
       else
-        write (*,'(a)')  '---- STEP LIMIT ----'
-      endif
-
-
-      write (*,'(/,2a,/)') ' Final energies ',
+        write (logString, '(a)')  '---- STEP LIMIT ----'
+      endif
+
+
+      write (logString, '(/,2a,/)') ' Final energies ',
      & '__________________________________________________'
 
@@ -197,12 +197,12 @@
       if (ireg.eq.0) then
 
-        write (*,'(a,e12.5,/,3(a,e11.4),/,2(a,e11.4))') ' Total: ',eysm,
-     &  '   Coulomb: ',eyel,' Lennard-Jones: ',eyvw,' HB: ',eyhb,
-     &  '   Variables: ',eyvr,'  Solvatation: ',eysl
+        write (logString, '(a,e12.5,/,3(a,e11.4),/,2(a,e11.4))') 
+     &  ' Total: ',eysm,'   Coulomb: ',eyel,' Lennard-Jones: ',eyvw,
+     &  ' HB: ',eyhb, '   Variables: ',eyvr,'  Solvatation: ',eysl
 
       else
 
-        write (*,'(a,e12.5,/,3(a,e11.4),/,3(a,e11.4))') ' Total: ',
-     &      wtey*eysm + wtrg*eyrg,
+        write (logString, '(a,e12.5,/,3(a,e11.4),/,3(a,e11.4))') 
+     &  ' Total: ', wtey*eysm + wtrg*eyrg,
      &  '   Coulomb: ',eyel,' Lennard-Jones: ',eyvw,' HB: ',eyhb,
      &  '   Variables: ',eyvr,'  Solvatation: ',eysl,
@@ -211,5 +211,5 @@
       endif
 
-      write (*,'(/,a,/)') ' Variables _________________'
+      write (logString, '(/,a,/)') ' Variables _________________'
 
       nv = 0
@@ -221,6 +221,6 @@
           if (abs(vr).gt.pi) vr=vr-sign(pi2,vr)
 
-          write (*,'(1x,a,1x,i4,f8.1,a,f5.1,a)') nmvr(i),nursvr(i),
-     &        vr*crd,'    (',abs(difang(vr,vlvro(i)))*crd,')'
+          write (logString, '(1x,a,1x,i4,f8.1,a,f5.1,a)') nmvr(i),
+     &        nursvr(i), vr*crd,'    (',abs(difang(vr,vlvro(i)))*crd,')'
 
           vlvr(i) = vr
@@ -231,13 +231,13 @@
       if (ireg.ne.0) then
 
-        write (*,'(/,a,/)') ' Global Variables ___________'
+        write (logString,'(/,a,/)') ' Global Variables ___________'
 
         do i=1,ntlml
-          write(*,*) ' Molecule #',i,' old          new'
+          write (logString, *) ' Molecule #',i,' old          new'
           do j=1,3
-            write(*,*) gbpro(j,i),' ',gbpr(j,i)
+            write (logString, *) gbpro(j,i),' ',gbpr(j,i)
           enddo
           do j=4,6
-            write(*,*) gbpro(j,i)*crd,' ',gbpr(j,i)*crd
+            write (logString, *) gbpro(j,i)*crd,' ',gbpr(j,i)*crd
           enddo
         enddo
@@ -245,13 +245,14 @@
       endif
 
-      write (*,'(/,2a)') ' Gradient ',
+      write (logString,'(/,2a)') ' Gradient ',
      & '______________________________________________________________'
 
-      write (*,'(8(1x,f8.3))') (gdvr(i),i=1,nv)
+      write (logString,'(8(1x,f8.3))') (gdvr(i),i=1,nv)
 
       if (ireg.ne.0) then
 
-        write (*,*) ' -------------- global variables ------------'
-        write (*,'(6(1x,f8.3))') (gdvr(i+nv),i=1,ngbvr)
+        write (logString, *) 
+     &             ' -------------- global variables ------------'
+        write (logString, '(6(1x,f8.3))') (gdvr(i+nv),i=1,ngbvr)
 
       endif
@@ -336,6 +337,6 @@
         esm=eysm
 
-        write (*,'(a,i5,a,2(e13.6,a))') ' Step ',nop,': energy ',esm
-     &                                 ,'  (',gdsmey,' )'
+        write (logString, '(a,i5,a,2(e13.6,a))') ' Step ',nop,
+     &                    ': energy ', esm, '  (',gdsmey,' )'
 
       else
@@ -355,6 +356,6 @@
         esm=wtey*eysm+wtrg*eyrg
 
-        write (*,'(a,i5,a,3(e13.6,a))') ' Step ',nop,': energy ',esm
-     &                                 ,'  (',gdsmey,',',gdsmrg,' )'
+        write (logString, '(a,i5,a,3(e13.6,a))') ' Step ',nop,
+     &             ': energy ',esm,'  (',gdsmey,',',gdsmrg,' )'
 
       endif
Index: minqsn.f
===================================================================
--- minqsn.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ minqsn.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -66,5 +66,6 @@
       integer ir
 
-
+      character(255) logString
+      
       parameter ( eps1=0.1d0,
      &            eps2=0.7d0,
@@ -142,5 +143,5 @@
 
     3 if (nfun.ge.maxfun) then
-!c        write (*,*) ' minfor> exceeded max. number of function calls'
+!c        write (logString, *) ' minfor> exceeded max. number of function calls'
         return
       endif
@@ -224,5 +225,6 @@
             goto 2
           else               ! rank of new matrix is deficient
-            write (*,*) ' minfor> rank of hessian < number of variables'
+            write (logString, *) 
+     &       ' minfor> rank of hessian < number of variables'
             return
           endif
Index: mklist.f
===================================================================
--- mklist.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ mklist.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -45,5 +45,5 @@
       ntlms=nmsml(nml)
       if (ntlms.eq.0) then
-        write (*,'(a,i4)')
+        write (logString, '(a,i4)')
      &           ' mklist> No mov. sets defined in molecule #',nml
         nvwml(nml)=0
@@ -107,5 +107,6 @@
                 n2i=n2i+1
                 if (n2i.gt.mx2) then
-                  write (*,'(a,i3,2a)')  ' mklist> Molecule # ',nml,
+                  write(logString, '(a,i3,2a)')  ' mklist> Molecule # ',
+     &                nml,
      &                ': too many atoms bound to ',nmat(iob)
                   stop
@@ -125,5 +126,6 @@
                     n2i=n2i+1
                     if (n2i.gt.mx2) then
-                      write (*,'(a,i3,a)')  ' mklist> Molecule # '
+                      write (logString, '(a,i3,a)')  
+     &                         ' mklist> Molecule # '
      &                         ,nml,': too many atoms in list L2I'
                       stop
@@ -141,5 +143,5 @@
                 n2i=n2i+1
                 if (n2i.gt.mx2) then
-                  write (*,'(a,i3,2a)')  ' mklist> Molecule # '
+                  write (logString, '(a,i3,2a)')  ' mklist> Molecule # '
      &             ,nml,': too many atoms bound to ',nmat(iob)
                   stop
@@ -165,5 +167,6 @@
               n1i=n1i+1
               if (n1i.gt.mxbd) then
-                write (*,'(a,i3,2a)')  ' mklist> Molecule # ',nml,
+                write (logString, '(a,i3,2a)')  ' mklist> Molecule # ',
+     &              nml,
      &             ': too many atoms bound to ',nmat(ib)
                 stop
@@ -177,6 +180,7 @@
                   n2i=n2i+1
                   if (n2i.gt.mx2) then
-                    write (*,'(a,i3,2a)')  ' mklist> Molecule # ',nml,
-     &              ': too many atoms bound to branches of ',nmat(ib)
+                    write (logString, '(a,i3,2a)')  
+     &                 ' mklist> Molecule # ',nml,
+     &                 ': too many atoms bound to branches of ',nmat(ib)
                     stop
                   endif
@@ -194,6 +198,7 @@
                 n2i=n2i+1
                 if (n2i.gt.mx2) then
-                  write (*,'(a,i3,a)')  ' mklist> Molecule # ',nml,
-     &            ': too many atoms n list L2I '
+                  write (logString, '(a,i3,a)')  ' mklist> Molecule # ',
+     &               nml,
+     &               ': too many atoms n list L2I '
                   stop
                 endif
@@ -207,5 +212,6 @@
 ! ________________________________________ Current 2ND list -> VdW-interact.
               if ((nvw+n2nd).gt.mxvw) then
-                write (*,'(a,i4,a,i5)') ' mklist> Molecule # ',nml,
+                write (logString, '(a,i4,a,i5)') ' mklist> Molecule # ',
+     &            nml,
      &            ': Number of vdw-domains > ',mxvw
                 stop
@@ -273,11 +279,11 @@
 ! ____________________________________________ Summary
 !      do i=ifiat,ilaat
-!        write (*,'(3a,i5,a)') ' ######## atom ',nmat(i),'(',i,')'
+!        write (logString, '(3a,i5,a)') ' ######## atom ',nmat(i),'(',i,')'
 !        iv1=ivwat1(i)
 !        iv2=ivwat2(i)
 !        if (iv1.le.iv2) then
-!          write(*,'(a)') ' ---> vdW :'
+!          write (logString, '(a)') ' ---> vdW :'
 !          do j=iv1,iv2
-!            write (*,'(i5,a,i5)') lvwat1(j),'-',lvwat2(j)
+!            write (logString, '(i5,a,i5)') lvwat1(j),'-',lvwat2(j)
 !          enddo
 !        endif
@@ -285,6 +291,6 @@
 !        i42=i14at2(i)
 !        if (i41.le.i42) then
-!          write(*,'(a)') ' ---> 1-4 :'
-!          write(*,'(10i5)') (l14at(j),j=i41,i42)
+!          write (logString, '(a)') ' ---> 1-4 :'
+!          write (logString, '(10i5)') (l14at(j),j=i41,i42)
 !        endif
 !      enddo
@@ -292,5 +298,5 @@
       return
 
-    1 write (*,'(a,i4,a,i5)') ' mklist> Molecule # ',nml,
+    1 write (logString, '(a,i4,a,i5)') ' mklist> Molecule # ',nml,
      &                     ': Number of 1-4 interactions > ',mx14
       stop
@@ -312,5 +318,5 @@
       implicit none
       integer mx, j, n, j1, l1, j2, l2, i1, i2, ja, k
-
+      character(255) logString
       dimension l1(mx),l2(mx)
 
@@ -334,5 +340,5 @@
               n=n+1
               if (n.gt.mx) then
-                write (*,'(a)') ' quench> too many sets'
+                write (logString, '(a)') ' quench> too many sets'
                 stop
               endif
Index: mulcan_par.f
===================================================================
--- mulcan_par.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ mulcan_par.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -83,6 +83,6 @@
 !
       eol = energy()
-      write (*,'(a,e12.5,/)')  'Energy of start configuration: ',eol
-      write(*,*)
+      write (logString, '(a,e12.5,/)')  'Energy of start configuration: ',eol
+      write (logString, *)
 
       call outpdb(1, 'start.pdb')
@@ -100,5 +100,5 @@
        muold = int(min(xmax,max(xmin,eol/ebin+sign(0.5d0,eol))))
        ihist(muold) = ihist(muold) + 1
-       write (*,*) nsw, eol, acz
+       write (logString, *) nsw, eol, acz
 !
 ! ITERATE MULTICANONICAL WEIGHTS EVERY NUP SWEEPS 
@@ -143,13 +143,13 @@
 ! FINAL OUTPUT:
       acz = acz/dble(nsw*nvr)
-      write(*,*) 'last energy',eol
-      write(*,*) 'aczeptance rate:',acz
+      write (logString, *) 'last energy',eol
+      write (logString, *) 'aczeptance rate:',acz
 !
 ! OUTPUT OF FINAL HISTOGRAM
 !
-      write(*,*) 'Histogram:'
+      write (logString, *) 'Histogram:'
       do i=kmin,kmax
        if(xhist(i).gt.0.0d0) then
-        write(*,*) i,xhist(i)
+        write (logString, *) i,xhist(i)
        end if
       end do
Index: mulcan_par_mod.f90
===================================================================
--- mulcan_par_mod.f90	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ mulcan_par_mod.f90	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -103,6 +103,6 @@
 
       eol = energy()
-      write (*,'(a,e12.5,/)')  'Energy of start configuration: ',eol
-      write(*,*)
+      write (logString, '(a,e12.5,/)')  'Energy of start configuration: ',eol
+      write (logString, *)
 
       call outpdb(1, 'start.pdb')
@@ -158,11 +158,11 @@
 ! Final output
       acz = acz/dble(nsw*nvr)
-      write(*,*) 'last energy',eol
-      write(*,*) 'acceptance rate:',acz
-
-      write(*,*) 'Histogram:'
+      write (logString, *) 'last energy',eol
+      write (logString, *) 'acceptance rate:',acz
+
+      write (logString, *) 'Histogram:'
       do i=kmin,kmax
        if(xhist(i).gt.0.0d0) then
-        write(*,*) i,xhist(i)
+        write (logString, *) i,xhist(i)
        end if
       end do
@@ -266,5 +266,5 @@
        end do
       end do
-      write(*,*) 'Number of contacts in reference conformation:',nci
+      write (logString, *) 'Number of contacts in reference conformation:',nci
 
 ! Read in fields with multicanonical parameter
@@ -282,5 +282,5 @@
         vlvr(iv) = x
        end do
-       write(*,*) 'Last iteration, energy:',nswm,eol_old
+       write (logString, *) 'Last iteration, energy:',nswm,eol_old
       else
 ! _________________________________ random start
@@ -294,13 +294,13 @@
 
       eol = energy()
-      write (*,'(e12.5,/)')  eol
+      write (logString, '(e12.5,/)')  eol
       call contacts(nhy,nhx,dham)
-      write(*,*) 'Number of contacts in start configuration:',nhy
-      write(*,*) 'Number of native contacts in start configuration:',   &
+      write (logString, *) 'Number of contacts in start configuration:',nhy
+      write (logString, *) 'Number of native contacts in start configuration:',   &
      &            nhx
       do i=1,nresi
-       write(*,'(62I1)') (ijcont(i,j), j=1,nresi)
-      end do
-      write(*,*)
+       write (logString, '(62I1)') (ijcont(i,j), j=1,nresi)
+      end do
+      write (logString, *)
 
 
@@ -355,11 +355,11 @@
 
       acz = acz/dble(nsw*nvr)
-      write(*,*) 'last energy',eol
-      write(*,*) 'acceptance rate:',acz
+      write (logString, *) 'last energy',eol
+      write (logString, *) 'acceptance rate:',acz
 
 ! WRITE DOWN (UN-REWEIGHTED) HISTOGRAM OF MULTICANONICAL SIMULATION
       do i=kmin,kmax
        if(xhist(i).gt.0.0d0) then
-        write(*,*) i,xhist(i)
+        write (logString, *) i,xhist(i)
        end if
       end do
Index: mulcan_sim.f
===================================================================
--- mulcan_sim.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ mulcan_sim.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -72,5 +72,5 @@
        end do
       end do
-      write(*,*) 'Number of contacts in reference conformation:',nci
+      write (logString, *) 'Number of contacts in reference conformation:',nci
 
 ! READ IN FIELDS WITH MULTICANONICAL PARAMETER
@@ -88,5 +88,5 @@
         vlvr(iv) = x
        end do
-       write(*,*) 'Last iteration, energy:',nswm,eol_old
+       write (logString, *) 'Last iteration, energy:',nswm,eol_old
       else
 ! _________________________________ random start
@@ -100,13 +100,13 @@
 !
       eol = energy()
-      write (*,'(e12.5,/)')  eol
+      write (logString, '(e12.5,/)')  eol
       call contacts(nhy,nhx,dham)
-      write(*,*) 'Number of contacts in start configuration:',nhy
-      write(*,*) 'Number of native contacts in start configuration:',
+      write (logString, *) 'Number of contacts in start configuration:',nhy
+      write (logString, *) 'Number of native contacts in start configuration:',
      &            nhx
       do i=1,nresi
-       write(*,'(62I1)') (ijcont(i,j), j=1,nresi)
-      end do
-      write(*,*)
+       write (logString, '(62I1)') (ijcont(i,j), j=1,nresi)
+      end do
+      write (logString, *)
 !
 
@@ -165,11 +165,11 @@
 ! FINAL OUTPUT:
       acz = acz/dble(nsw*nvr)
-      write(*,*) 'last energy',eol
-      write(*,*) 'aczeptance rate:',acz
+      write (logString, *) 'last energy',eol
+      write (logString, *) 'aczeptance rate:',acz
 
 ! WRITE DOWN (UN-REWEIGHTED) HISTOGRAM OF MULTICANONICAL SIMULATION
       do i=kmin,kmax
        if(xhist(i).gt.0.0d0) then
-        write(*,*) i,xhist(i)
+        write (logString, *) i,xhist(i)
        end if
       end do
Index: nursvr.f
===================================================================
--- nursvr.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ nursvr.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -33,5 +33,6 @@
       enddo
 
-      write (*,'(a,i5)') ' nursvr > Cannot find variable # ',ivr
+      write (logString, '(a,i5)') ' nursvr > Cannot find variable # '
+     &   ,ivr
       stop
 
@@ -70,5 +71,5 @@
       enddo
 
-      write (*,'(a,i5)') ' nursat > Cannot find atom # ',iat
+      write (logString, '(a,i5)') ' nursat > Cannot find atom # ',iat
       stop
 
Index: opeflx.f
===================================================================
--- opeflx.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ opeflx.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -60,5 +60,5 @@
       ntlvr=nvrml(nml)
       if (ntlvr.eq.0) then
-        write (*,'(a,i4)')
+        write (logString, '(a,i4)')
      &           ' opeflx> No variables defined in molecule #',nml
         return
@@ -504,5 +504,5 @@
       gda=gdeyvr(iv)         ! analytical der.
 
-      write (*,'(1x,2a,2(e12.6,a))') nmvr(iv),': ',gda,' (',
+      write (logString, '(1x,2a,2(e12.6,a))') nmvr(iv),': ',gda,' (',
      &       abs(gda-gdn),')'
 
Index: opereg.f
===================================================================
--- opereg.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ opereg.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -54,5 +54,5 @@
       ntlvr=nvrml(nml)
       if (ntlvr.eq.0) then
-        write (*,'(a,i4)')
+        write (logString, '(a,i4)')
      &           ' opereg> No variables defined in molecule #',nml
         return
@@ -330,5 +330,5 @@
         gdn = ( enyreg(nml) - eyrg ) / del
 
-        write (*,*) ' Gb. var #',(ii+i),': ',gdeygb(ii+i),gdn,
+        write (logString, *) ' Gb. var #',(ii+i),': ',gdeygb(ii+i),gdn,
      &                                   abs(gdn-gdeygb(ii+i))
 ! ----------------------------- restore
@@ -382,5 +382,5 @@
       gda=gdeyrg(iv)         ! analytical der.
 
-      write (*,'(1x,2a,2(e12.6,a))') nmvr(iv),': ',gda,' (',
+      write (logString, '(1x,2a,2(e12.6,a))') nmvr(iv),': ',gda,' (',
      &       abs(gda-gdn),')'
 
Index: opeshe.f
===================================================================
--- opeshe.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ opeshe.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -51,5 +51,5 @@
       ntlvr=nvrml(nml)
       if (ntlvr.eq.0) then
-        write (*,'(a,i4)')
+        write (logString, '(a,i4)')
      &           ' opeshe> No variables defined in molecule #',nml
         return
@@ -377,5 +377,5 @@
       gda=gdeyvr(iv)         ! analytical der.
 
-      write (*,'(1x,2a,2(e12.6,a))') nmvr(iv),': ',gda,' (',
+      write (logString, '(1x,2a,2(e12.6,a))') nmvr(iv),': ',gda,' (',
      &       abs(gda-gdn),')'
 
Index: opesol.f
===================================================================
--- opesol.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ opesol.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -51,5 +51,5 @@
       ntlvr=nvrml(nml)
       if (ntlvr.eq.0) then
-        write (*,'(a,i4)')
+        write (logString, '(a,i4)')
      &           ' opesol> No variables defined in molecule #',nml
         return
@@ -283,5 +283,5 @@
       gdn=(eynw-eysl)/del    ! numerical derivative
 
-      write (*,'(1x,2a,2(e12.6,a))') nmvr(iv),': ',gda,' (',
+      write (logString, '(1x,2a,2(e12.6,a))') nmvr(iv),': ',gda,' (',
      &       abs(gda-gdn),')'
 
Index: outpdb.f
===================================================================
--- outpdb.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ outpdb.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -43,5 +43,5 @@
 
       if (nml.lt.0.or.nml.gt.ntlml) then
-        write(*,*) ' outpdb>  No such molecule #',nml,' !'
+        write (logString, *) ' outpdb>  No such molecule #',nml,' !'
         return
       elseif (nml.gt.0) then
Index: outvar.f
===================================================================
--- outvar.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ outvar.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -30,5 +30,5 @@
 
       if (nml.lt.0.or.nml.gt.ntlml) then
-        write(*,*) ' outvar>  No such molecule #',nml,' !'
+        write (logString, *) ' outvar>  No such molecule #',nml,' !'
         return
       endif
@@ -57,5 +57,5 @@
 
       if (nml.lt.0.or.nml.gt.ntlml) then
-        write(*,*) ' outvbs>  No such molecule #',nml,' !'
+        write (logString, *) ' outvbs>  No such molecule #',nml,' !'
         return
       elseif (nml.gt.0) then
Index: partem_p.f
===================================================================
--- partem_p.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ partem_p.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -95,6 +95,6 @@
       trackID = 1
       odd = 1
-      write (*,*) 'Starting parallel tempering.'
-      write (*,*) 'parameters, ',switch,newsta,nmes,nswp,nmes,
+      write (logString, *) 'Starting parallel tempering.'
+      write (logString, *) 'parameters, ',switch,newsta,nmes,nswp,nmes,
      &            rep_id, num_rep, partem_comm, myrank
       call flush(6)
@@ -167,5 +167,5 @@
             do i=1,num_rep
                read(13,*) j,inode(i),intem(i),yol(i),e_minp(i),h_maxp(i)
-               write (*,*) "par_R.in:",i,j
+               write (logString, *) "par_R.in:",i,j
             end do
             jold=(iold/nmes)*num_rep
@@ -178,5 +178,5 @@
      &              ,idum1, idum2, idum3, e_min
      &              ,dummy, dummy
-               write (*,*) i
+               write (logString, *) i
                call flush(6)
             end do
@@ -197,12 +197,12 @@
       e_min = e_minp(rep_id+1)
       h_max = h_maxp(rep_id+1)
-      write (*,*) "E_min=",e_min," for ", rep_id + 1 
+      write (logString, *) "E_min=",e_min," for ", rep_id + 1 
       eol=energy()
       if(.not.newsta.and.abs(yol(rep_id + 1) - eol).gt.0.1) then
-         write(*,*) rep_id, ' Warning: yol(rep_id).ne.eol:'
-         write(*,*) rep_id, yol(rep_id + 1), eol
+         write (logString, *) rep_id, ' Warning: yol(rep_id).ne.eol:'
+         write (logString, *) rep_id, yol(rep_id + 1), eol
       endif
 !     Start of simulation
-      write (*,*) '[',rep_id, myrank, beta, partem_comm,
+      write (logString, *) '[',rep_id, myrank, beta, partem_comm,
      &            '] Energy before equilibration:', eol
 !     =====================Equilibration by canonical Metropolis
@@ -211,5 +211,6 @@
       end do
       CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
-      write (*,*) '[',rep_id,'] Energy after equilibration:', eol
+      write (logString, *) '[',rep_id,'] Energy after equilibration:', 
+     &   eol
       call flush(6)
 !     
@@ -222,5 +223,5 @@
          eol0 = eol
          if (myrank.eq.0.and.rep_id.eq.0) then
-            write (*,*) "Finished sweep", nsw
+            write (logString, *) "Finished sweep", nsw
             call flush(6)
          endif
@@ -494,9 +495,9 @@
      &                IERR)
                if (myrank.ne.0) then
-!                  write (*,*) '[', myrank,'] Missed', randomCount, 
+!                  write (logString, *) '[', myrank,'] Missed', randomCount, 
 !     &                            'random numbers.'
                   do i = 1, randomCount
                      rd = grnd()
-!                     write (*,*) '[', myrank,'] rd=', rd
+!                     write (logString, *) '[', myrank,'] rd=', rd
                   enddo
                endif
@@ -519,5 +520,5 @@
       e_final=energy()
       if (partem_comm.ne.MPI_COMM_NULL) then
-         write (*,*) rep_id, ' E_final', e_final
+         write (logString, *) rep_id, ' E_final', e_final
       endif
       eol0 = eol
Index: partem_s.f
===================================================================
--- partem_s.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ partem_s.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -148,6 +148,6 @@
                CALL METROPOLIS(energ,acz,can_weight)
             end do
-            write(*,*) 'Start energy after equilibration for replica:',
-     &                 k, energ
+            write (logString, *) 
+     &         'Start energy after equilibration for replica:',k, energ
             do i=1,nvr
                iv = idvr(i)
Index: pdbread.f
===================================================================
--- pdbread.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ pdbread.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -32,5 +32,5 @@
       dimension cor(3)
       character atm*4,rsn*3,rsno*3,chn,chno,
-     &          rsid*5,rsido*5,line*132
+     &          rsid*5,rsido*5,line*132, logString*255
 
       natp=0
@@ -48,5 +48,5 @@
         lunpdb = 99
       else
-        write (*,'(a)')
+        write (logString, '(a)')
      &    ' pdbread> empty file name to read pdb-structure'
 
@@ -57,5 +57,5 @@
 
       if (io.le.0) then
-        write (*,'(a,/,a)')
+        write (logString, '(a,/,a)')
      &  ' pdbread> ERROR opening file to read pdb-structure: ',
      &  pdbfil(1:iendst(pdbfil))
@@ -70,5 +70,5 @@
 
       if ( line(17:17).ne.' ' )  then
-        write (*,'(a,/,a,/,a,/,2a)')
+        write (logString, '(a,/,a,/,a,/,2a)')
      &  ' pdbread> found alternate atom location: ',
      &  '                !',
@@ -86,5 +86,5 @@
 
       if ((natp+1).gt.MXATP) then
-        write (*,'(a,i5,a,/,a)')
+        write (logString, '(a,i5,a,/,a)')
      &  ' pdbread>  >MXATP (',MXATP,') ATOM lines in PDB file ',
      &  pdbfil(1:iendst(pdbfil))
@@ -97,5 +97,5 @@
 
         if ((nchp+1).gt.MXCHP) then
-          write (*,'(a,i3,a,/,a)')
+          write (logString, '(a,i3,a,/,a)')
      &    ' pdbread>  >MXCHP (',MXCHP,') chains in PDB file ',
      &    pdbfil(1:iendst(pdbfil))
@@ -106,5 +106,5 @@
 
         if ((nrsp+1).gt.MXRSP) then
-          write (*,'(a,i3,a,/,a)')
+          write (logString, '(a,i3,a,/,a)')
      &    ' pdbread>  >MXRSP (',MXRSP,') residues in PDB file ',
      &    pdbfil(1:iendst(pdbfil))
@@ -141,5 +141,5 @@
 
         if ((nrsp+1).gt.MXRSP) then
-          write (*,'(a,i3,a,/,a)')
+          write (logString, '(a,i3,a,/,a)')
      &    ' pdbread>  >MXRSP (',MXRSP,') residues in PDB file ',
      &    pdbfil(1:iendst(pdbfil))
@@ -172,5 +172,5 @@
       goto 1
 
-    2 write (*,'(a,/,a,/,2a)')
+    2 write (logString, '(a,/,a,/,2a)')
      &  ' pdbread> ERROR reading ATOM line ',
      &  line(:l),
@@ -195,5 +195,5 @@
       else
 
-        write (*,'(a,/,a)')
+        write (logString, '(a,/,a)')
      &  ' pdbread> NO atom coordinates selected from file ',
      &  pdbfil(1:iendst(pdbfil))
@@ -245,5 +245,5 @@
         nml=nml+1
         if (nml.gt.mxml) then
-          write(*,'(a,i4,2a)')' pdbvars> NUMBER of chains > '
+          write (logString, '(a,i4,2a)')' pdbvars> NUMBER of chains > '
      &                          ,mxml,' in ',' ?'
           stop
@@ -266,6 +266,6 @@
 
           if (nrs.gt.mxrs) then
-            write(*,'(a,i4,2a)') ' pdbvars> NUMBER of residues > '
-     &                       ,mxrs,' in ',' ?'
+            write (logString, '(a,i4,2a)') 
+     &         ' pdbvars> NUMBER of residues > ', mxrs, ' in ',' ?'
             stop
           endif
@@ -421,6 +421,6 @@
             zat(ii)=zatp(ix)
           else
-            write(*,'(3a)') ' pdbvars> missing PDB atom ',nmat(ii),
-     &       ' is ref. point for SMMP - cannot proceed !'
+            write (logString, '(3a)') ' pdbvars> missing PDB atom ',
+     &         nmat(ii), ' is ref. point for SMMP - cannot proceed !'
           endif
         enddo
@@ -436,6 +436,6 @@
 ! ++++++++++
 
-       write(*,*) ' '
-       write(*,*) ' Initial RMSD ',rmsd
+       write (logString, *) ' '
+       write (logString, *) ' Initial RMSD ',rmsd
 
       enddo ! chains(molecules)
@@ -485,5 +485,5 @@
             enddo
 
-!            write(*,'(8a)') ' pdbvars> ',atm,' not found in '
+!            write (logString, '(8a)') ' pdbvars> ',atm,' not found in '
 !     #       ,chnp(nc),' ',rsidp(irs),' ',rsnmp(irs)
 
Index: redseq.f
===================================================================
--- redseq.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ redseq.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -36,5 +36,5 @@
 
       if (iopfil(lunseq,seqfil,'old','formatted').le.izero) then
-        write (*,'(a,/,a,i3,2a)') 
+        write (logString, '(a,/,a,i3,2a)') 
      &    ' redseq> ERROR opening sequence file:',
      &      ' LUN=',lunseq,' FILE=',seqfil(1:iendst(seqfil))
@@ -67,5 +67,5 @@
 
             if ((nrs-irsml1(ntlml)+1).eq.0) then
-              write(*,'(2a)') ' redseq> IGNORE molecule: ',
+              write (logString, '(2a)') ' redseq> IGNORE molecule: ',
      &                        nmml(ntlml)(1:iendst(nmml(ntlml)))
               ntlml=ntlml-1
@@ -74,5 +74,6 @@
           ntlml=ntlml+1
           if (ntlml.gt.mxml) then
-            write(*,'(a,i4,2a)')' redseq> NUMBER of molecules > '
+            write (logString, '(a,i4,2a)') 
+     &                          ' redseq> NUMBER of molecules > '
      &                          ,mxml,' in ',seqfil(1:iendst(seqfil))
             close(lunseq)
@@ -124,5 +125,5 @@
 
           if (id.gt.4) then
-            write (*,'(4a)') ' redseq> INVALID residue NAME >',
+            write (logString, '(4a)') ' redseq> INVALID residue NAME >',
      &                       line(ib:ii),'< in ',
      &      seqfil(1:iendst(seqfil))
@@ -133,5 +134,6 @@
             nrs=nrs+1
             if (nrs.gt.mxrs) then
-              write(*,'(a,i4,2a)') ' redseq> NUMBER of residues > '
+              write (logString, '(a,i4,2a)') 
+     &                       ' redseq> NUMBER of residues > '
      &                       ,mxrs,' in ',seqfil(1:iendst(seqfil))
               close(lunseq)
@@ -162,5 +164,5 @@
 
       if (nrs.eq.0) then
-        write (*,'(2a)') ' redseq> no residues found in ',
+        write (logString, '(2a)') ' redseq> no residues found in ',
      &                   seqfil(1:iendst(seqfil))
         stop
@@ -179,10 +181,11 @@
           if (i.eq.ntlml) then   ! Check last molecule
             if ((nrs-ifirs+1).eq.0) then
-              write(*,'(2a)') ' redseq> IGNORE molecule '
+              write (logString, '(2a)') ' redseq> IGNORE molecule '
      &                        ,nmml(ntlml)(1:iendst(nmml(ntlml)))
               ntlml=ntlml-1
               if (ntlml.eq.0) then
-                write (*,'(2a)') ' redseq> no residues found in ',
-     &          seqfil(1:iendst(seqfil))
+                write (logString, '(2a)') 
+     &             ' redseq> no residues found in ',
+     &             seqfil(1:iendst(seqfil))
                 stop
               endif
@@ -192,8 +195,8 @@
           endif
 
-!          write (*,'(/,a,i4,2a)') ' redseq> ',irsml2(i)-irsml1(i)+1,
+!          write (logString, '(/,a,i4,2a)') ' redseq> ',irsml2(i)-irsml1(i)+1,
 !     &           ' residue(s) in molecule: ',
 !     &           nmml(i)(1:iendst(nmml(i)))
-!          write (*,'(15(1x,a))') (seq(j),j=irsml1(i),irsml2(i))
+!          write (logString, '(15(1x,a))') (seq(j),j=irsml1(i),irsml2(i))
 
         enddo
@@ -203,5 +206,5 @@
 ! _______________________________________________ error
 
-    4 write (*,'(a,i4,2a)') ' redseq> ERROR reading line No. ',nln,
+    4 write (logString, '(a,i4,2a)') ' redseq> ERROR reading line No. ',nln,
      & ' in ',seqfil(1:iendst(seqfil))
       close(lunseq)
Index: redstr.f
===================================================================
--- redstr.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ redstr.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -30,4 +30,5 @@
 
       character spr,blnk,str*(*),strn*(*)
+      character*255 logString
       
       data blnk/' '/
@@ -36,5 +37,5 @@
 
       if (spr.eq.blnk) then
-        write (*,*) ' extstr> Separator should not be blank'
+        write (logString, *) ' extstr> Separator should not be blank'
         stop
       endif
@@ -79,5 +80,6 @@
       return
 ! ______________________________________________________________ Error
-    1 write (*,'(a)') ' extstr> Substring to be extracted is too long !'
+    1 write (logString, '(a)') 
+     &   ' extstr> Substring to be extracted is too long !'
       stop
 
Index: redvar.f
===================================================================
--- redvar.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ redvar.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -101,5 +101,5 @@
       ntlvr=ivrml1(ntlml)+nvrml(ntlml)-1
       if (ntlvr.eq.0) then
-        write (*,*) ' redvar> No variables defined in molecule(s)'
+        write(logString,*)'redvar> No variables defined in molecule(s)'
         return
       endif
@@ -108,5 +108,5 @@
       io=iopfil(lunvar,varfil,'old','formatted')
       if (io.eq.0) then
-        write (*,'(a,/,a,i3,2a)') 
+        write (logString, '(a,/,a,i3,2a)') 
      &    ' redvar> ERROR opening file to set variables:',
      &    ' LUN=',lunvar,' FILE=',varfil(1:iendst(varfil))
@@ -136,5 +136,5 @@
           vr=blat(iatvr(i))
         else
-          write(*,*) 'redvar>  unknown variable type: ',it,' !'
+          write(logString,*) 'redvar>  unknown variable type: ',it,' !'
           stop
         endif
@@ -159,5 +159,6 @@
             if (iredin(lincmd,nml).le.0.or.
      &          nml.le.0.or.nml.gt.ntlml) then
-              write (*,*) 'redvar> ','Incorrect molecule number >',
+              write (logString, *) 
+     &                 'redvar> ','Incorrect molecule number >',
      &                    lincmd(1:l),'<  Must be in range [1,',
      &                    ntlml,'] !'
@@ -310,5 +311,5 @@
 
                         if (inum.le.0.or.inum.gt.ntlml) then
-                          write (*,*) ' # 1: ',inum
+                          write (logString, *) ' # 1: ',inum
                           goto 104
                         endif
@@ -324,5 +325,5 @@
 
                             if (k.lt.nfi.or.k.gt.irsml2(j)) then
-                              write (*,*) ' # 2: ',inum
+                              write (logString, *) ' # 2: ',inum
                               goto 104
                             endif 
@@ -343,5 +344,5 @@
                                 if (l.lt.nfi.or.
      &                              l.gt.nfi+nvrrs(k)-1) then
-                                  write (*,*) ' # 3: ',inum
+                                  write (logString, *) ' # 3: ',inum
                                   goto 104
                                 endif
@@ -427,5 +428,5 @@
 
                     if (ide.eq.1.and.ihy.eq.ide) then
-                      write (*,*) ' # 4: ',ide
+                      write (logString, *) ' # 4: ',ide
                       goto 104
                     endif
@@ -441,5 +442,5 @@
                       if (iredin(linh,ibz).le.0.or.ibz.le.0)
      &                  then
-                        write (*,*) ' # 5 '
+                        write (logString, *) ' # 5 '
                         goto 104
                       endif
@@ -456,5 +457,5 @@
                       if (iredin(linh,iez).le.0.or.iez.le.0.or.
      &                  iez.lt.ibz) then
-                        write (*,*) ' # 6 '
+                        write (logString, *) ' # 6 '
                         goto 104
                       endif
@@ -524,5 +525,5 @@
                 endif
               enddo
-              if (.not.did) write (*,'(3a)')
+              if (.not.did) write (logString, '(3a)')
      &          ' redvar> No variables affected by command >',
      &          lincmd(1:ice),'<'
@@ -531,5 +532,5 @@
               ll1=ibegst(linfld(nfld))
               ll2=iendst(linfld(nfld))
-              write (*,*) 'll1,ll2, linfld(nfld): ',ll1,ll2,
+              write (logString, *) 'll1,ll2, linfld(nfld): ',ll1,ll2,
      &           '>',linfld(nfld)(ll1:ll2),'<'
 
@@ -559,5 +560,6 @@
              enddo
 
-             write (*,'(3a,/,1x,5(a,2x),a)') ' redvar> ',nmml(i)(1:ie),
+             write (logString, '(3a,/,1x,5(a,2x),a)') ' redvar> ',
+     &                               nmml(i)(1:ie),
      &                                    ' with global parameters:',
      &                              (strg(k)(ibegst(strg(k)):),k=1,6)
@@ -579,9 +581,9 @@
                 did=.true.
                 if (fxvr(iv)) then
-                  write (*,'(3a,i4,1x,4a,f10.3,a)') ' redvar> ',
+                  write(logString,'(3a,i4,1x,4a,f10.3,a)') ' redvar> ',
      &                nmml(i)(1:ie),': residue ',j-jb,seq(j),
      &                ': ',nmvr(iv),' set ',vlvrx(iv),'   Fixed'
                 else
-                  write (*,'(3a,i4,1x,4a,f10.3)') ' redvar> ',
+                  write (logString, '(3a,i4,1x,4a,f10.3)') ' redvar> ',
      &                nmml(i)(1:ie),': residue ',j-jb,seq(j),
      &                ': ',nmvr(iv),' set ',vlvrx(iv)
@@ -597,10 +599,10 @@
           enddo
           if (did) then
-            if (in.gt.0) write (*,'(3a,i5,a)')
+            if (in.gt.0) write (logString, '(3a,i5,a)')
      &        ' redvar> Molecule ',nmml(i)(1:ie),': ',in,
      &        ' variable(s) remain unchanged'
             call setvar(iml,vlvrx)
           else
-            write (*,'(3a)') ' redvar> Molecule ',
+            write (logString, '(3a)') ' redvar> Molecule ',
      &        nmml(i)(1:ie),': No internal variables changed'
           endif
@@ -610,24 +612,26 @@
       return
 ! ____________________________________________________________ Errors
-  100 write (*,'(3a)') ' redvar> Cannot interpret command >',
+  100 write (logString, '(3a)') ' redvar> Cannot interpret command >',
      &                 lincmd(1:ice),'<'
       close(lunvar)
       stop
-  101 write (*,'(a,i5,a)') ' redvar> Command number ',ncmd,' reached'
+  101 write (logString, '(a,i5,a)') ' redvar> Command number ',ncmd,
+     &       ' reached'
       close(lunvar)
       stop
-  102 write (*,'(3a)') ' redvar> Cannot read value from >',
+  102 write (logString, '(3a)') ' redvar> Cannot read value from >',
      &                   lincmd(1:ice),'<'
       close(lunvar)
       stop
-  103 write (*,'(a,i3,3a)') ' redvar> Cannot read >',maxide,
+  103 write (logString, '(a,i3,3a)') ' redvar> Cannot read >',maxide,
      &         ' identifiers from >',linfld(ii)(1:ife),'<'
       close(lunvar)
       stop
-  104 write (*,'(5a)') ' redvar> Error in identifier >',
+  104 write (logString, '(5a)') ' redvar> Error in identifier >',
      &            linide(1:ide),'< of command >',lincmd(1:ice),'<'
       close(lunvar)
       stop
-  105 write (*,'(a,/,a,/,2a,/)') ' redvar> line with global paramters:',
+  105 write (logString, '(a,/,a,/,2a,/)') 
+     &                           ' redvar> line with global paramters:',
      &                           line(1:ile),' must contain 6 floating',
      &                           ' point numbers separated by commas !'
@@ -635,5 +639,6 @@
       stop
 
-  106 write (*,'(a,/,a,/,2a,/)') ' redvar> line with global paramters:',
+  106 write (logString, '(a,/,a,/,2a,/)') 
+     &                           ' redvar> line with global paramters:',
      &                           line(1:ile),' angles must be inside ',
      &'ranges [-180,180], [-90,90], and [-180,180] Deg., respectively !'
Index: regul.f
===================================================================
--- regul.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ regul.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -44,5 +44,5 @@
       wtey = 0.d0
 
-      write(*,'(/,a,2(a,f4.2),/)')
+      write (logString, '(/,a,2(a,f4.2),/)')
      &  ' ====================== Regularization only',
      &  '   Wt(energy) = ',wtey,'  Wt(regul.) = ',wtrg
@@ -50,14 +50,14 @@
       call minim(2, nsteps, acc)
 
-      write(*,*) ' '
-      write(*,*) ' ---------- contacts after 1st regularization'
-      write(*,*) ' '
+      write (logString, *) ' '
+      write (logString, *) ' -------- contacts after 1st regularization'
+      write (logString, *) ' '
       call cnteny(nml)
-      write(*,*) ' '
+      write (logString, *) ' '
 
       nrs = irsml2(nml)-irsml1(nml)+1
       call rmsdopt(nml,1,nrs,ixatp,xatp,yatp,zatp,0,rm,av1,av2,rmsd)
 
-      write(*,*) ' RMSD = ',rmsd
+      write (logString, *) ' RMSD = ',rmsd
 
 ! --------------------------------------- fix vars. defined in PDB
@@ -70,5 +70,5 @@
       ireg = 0
 
-      write(*,'(/,a,2(a,f4.2),/)')
+      write (logString, '(/,a,2(a,f4.2),/)')
      &  ' ====================== Internal Energy for Hydrogens only',
      &  '   Wt(energy) = ',wtey,'  Wt(regul.) = ',wtrg
@@ -76,7 +76,7 @@
       call minim(1, nsteps, acc)
 
-      write(*,*) ' '
-      write(*,*) ' ---------- contacs after Emin. for hydrogens'
-      write(*,*) ' '
+      write (logString, *) ' '
+      write (logString, *) ' -------- contacs after Emin. for hydrogens'
+      write (logString, *) ' '
       call cnteny(nml)
 
@@ -98,5 +98,5 @@
         wtey = 1.d0 - wtrg
 
-        write(*,'(/,a,i2,2(a,e11.3),/)')
+        write (logString, '(/,a,i2,2(a,e11.3),/)')
      &    ' ================ Minimization #',it,
      &        '   Wt(energy) = ',wtey,'  Wt(regul.) = ',wtrg
@@ -107,12 +107,12 @@
         call rmsdopt(nml,1,nrs,ixatp,xatp,yatp,zatp,0,rm,av1,av2,rmsd)
 
-        write(*,*) ' '
-        write(*,*) ' RMSD = ',rmsd
+        write (logString, *) ' '
+        write (logString, *) ' RMSD = ',rmsd
 
       enddo
 
-      write(*,*) ' '
-      write(*,*) ' ---------- contacts after full regularization'
-      write(*,*) ' '
+      write (logString, *) ' '
+      write (logString, *) ' ------- contacts after full regularization'
+      write (logString, *) ' '
       call cnteny(nml)
 
@@ -120,5 +120,6 @@
 
 ! Output of dihedral angles of the regularized structure
-      write(*,*) 'Dihedral angles of the regularized structure;'
+      write (logString, *) 
+     &   'Dihedral angles of the regularized structure;'
       call outvar(nml, 'regd.var')
 
Index: rgyr.f
===================================================================
--- rgyr.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ rgyr.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -48,5 +48,5 @@
 
       if (nat.le.0) then
-        write (*,'(a,i4)')
+        write (logString, '(a,i4)')
      &     ' rgyr> No atoms found for molecule #',nml
         return 
Index: rmsdfun.f
===================================================================
--- rmsdfun.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ rmsdfun.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -82,5 +82,6 @@
 
       if (nml.lt.1.or.nml.gt.ntlml) then
-        write(*,*) ' rmsdopt>  Sorry, there is no molecule #',nml
+        write (logString, *) ' rmsdopt>  Sorry, there is no molecule #',
+     &     nml
         stop
       endif
@@ -137,5 +138,5 @@
 
       if (n.lt.3) then
-        write(*,*) ' rmsdopt>  <3 atoms selected !'
+        write (logString, *) ' rmsdopt>  <3 atoms selected !'
         stop
       endif
@@ -290,5 +291,5 @@
       integer n,nrot,i,ip,iq,j
 
-
+      character(255) logString
       real*8 a(n,n),d(n),v(n,n),
      &       c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX),smeps
@@ -403,5 +404,5 @@
       enddo
 
-      write(*,*) ' jacobi> too many iterations'
+      write (logString, *) ' jacobi> too many iterations'
       stop
 
Index: setmvs.f
===================================================================
--- setmvs.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ setmvs.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -77,5 +77,5 @@
 
       if (ntlvr.eq.0) then
-        write (*,'(a,i4)')
+        write (logString, '(a,i4)')
      &           ' setmvs> No variables defined in molecule #',nml
         nmsml(nml)=0
@@ -176,5 +176,5 @@
               if (i1.ne.0) then  ! combine ranges
                 if (j.gt.(i2+1).or.k.lt.(i1-1)) then
-                  write (*,'(3a,/,2a,i4,a,i3)') 
+                  write (logString, '(3a,/,2a,i4,a,i3)') 
      &             ' setmvs> Cannot combine disjunct ranges of atom',
      &             ' indices for torsion ',nmvr(iv),' in residue ',
@@ -198,5 +198,5 @@
 
         if ((nms+1).gt.mxms) then
-          write (*,'(a,i4,a,i5)') ' setmvs> Molecule # ',nml,
+          write (logString, '(a,i4,a,i5)') ' setmvs> Molecule # ',nml,
      &    ': Number of moving sets > ',mxms
           stop
@@ -231,5 +231,6 @@
                   nms=nms+1
                   if (nms.gt.mxms) then
-                    write (*,'(a,i4,a,i5)') ' setmvs> Molecule # ',
+                    write (logString, '(a,i4,a,i5)') 
+     &                ' setmvs> Molecule # ',
      &               nml,': Number of moving sets > ',mxms
                      stop
@@ -320,5 +321,6 @@
             nad=nad+1
             if (nad.gt.mxvr) then
-              write (*,'(a,i4,a,i5)') ' setmvs> Molecule # ',nml,
+              write (logString, '(a,i4,a,i5)') ' setmvs> Molecule # ',
+     &                         nml,
      &                         ': Number of added variables > ',mxvr
               stop
@@ -335,5 +337,6 @@
                 nad=nad+1
                 if (nad.gt.mxvr) then
-                  write (*,'(a,i4,a,i5)') ' setmvs> Molecule # ',nml,
+                  write (logString, '(a,i4,a,i5)') 
+     &                         ' setmvs> Molecule # ',nml,
      &                         ': Number of added variables > ',mxvr
                   stop
@@ -349,5 +352,6 @@
               nad=nad+1
               if (nad.gt.mxvr) then
-                write (*,'(a,i4,a,i5)') ' setmvs> Molecule # ',nml,
+                write (logString, '(a,i4,a,i5)') ' setmvs> Molecule # ',
+     &                       nml,
      &                       ': Number of added variables > ',mxvr
                 stop
@@ -375,14 +379,14 @@
 !            i2=latms2(i)
 !            if (i.eq.i1s) then
-!              write (*,'(a,i3,7a,i4,3a,i4,a)') 'res # ',nursvr(iv),
+!              write (logString, '(a,i3,7a,i4,3a,i4,a)') 'res # ',nursvr(iv),
 !     #        ' var: ',nmvr(iv),' base:',nmat(ib),'    atoms= ',
 !     #        nmat(i1),'(',i1,') - ',nmat(i2),'(',i2,')'
 !            else
-!              write (*,'(39x,2a,i4,3a,i4,a)')
+!              write (logString, '(39x,2a,i4,3a,i4,a)')
 !     #        nmat(i1),'(',i1,') - ',nmat(i2),'(',i2,')'
 !            endif
 !          enddo
 !        else
-!          write (*,'(a,i3,5a)') 'res # ',nursvr(iv),
+!          write (logString, '(a,i3,5a)') 'res # ',nursvr(iv),
 !     #    ' var: ',nmvr(iv),' base:',nmat(ib),'  No atoms'
 !        endif
@@ -390,8 +394,8 @@
 !        i2a=iadvr2(iv)
 !        if (i1a.le.i2a) then
-!          write (*,'(a,30(1x,a))') ' Depending variables:',
+!          write (logString, '(a,30(1x,a))') ' Depending variables:',
 !     #                    (nmvr(ladvr(i)),i=i1a,i2a)
 !        else
-!          write (*,'(a)') ' No dep. variables'
+!          write (logString, '(a)') ' No dep. variables'
 !        endif
 !      enddo
@@ -400,5 +404,5 @@
       return
 
-    6 write (*,'(a,i4,/,2(a,i5),a)') 
+    6 write (logString, '(a,i4,/,2(a,i5),a)') 
      & ' setmvs> Error in atom numbering of molecule # ',nml,
      & ': atom ranges for variables # ',iv,' and # ',jv,
Index: utilities.f
===================================================================
--- utilities.f	(revision 32289cdc085b7eebe6f6df5975767c27f3bf61c7)
+++ utilities.f	(revision 38d77eb4e27a6b22e4388b28629620f02aec9598)
@@ -66,6 +66,7 @@
       end do
       ipps = isum / num_ppr
-      write (*,*) "Total number of interactions:", isum
-      write (*,*) "Average # of interactions per processor", ipps
+      write (logString, *) "Total number of interactions:", isum
+      write (logString, *) "Average # of interactions per processor", 
+     &   ipps
 
       totalct = 0
@@ -101,8 +102,8 @@
      &            .lt.abs(totalct + isum - itarget))) then
                workPerProcessor(nml, irank) = io + 1
-!                write (*,*) io + 1, totalct, itarget
+!                write (logString, *) io + 1, totalct, itarget
             else
                workPerProcessor(nml, irank) = io
-!                write (*,*) io, totalct + isum, itarget
+!                write (logString, *) io, totalct + isum, itarget
             end if
             irank = irank + 1
@@ -127,5 +128,5 @@
 !     rank = 11
 !     fileName = fileNameMP('base_0000.dat', 6, 9, rank)
-!     write (*,*), fileName
+!     write (logString, *), fileName
 !     \endcode
 !     will output base_0011.dat.
