Fixed eello5, eello6, eturn6, and shortrange RESPA

diff --git a/bin/unres/MD/unres-mult-symetr_ifort_MPICH_GAB.exe b/bin/unres/MD/unres-mult-symetr_ifort_MPICH_GAB.exe
index 931588b..9d96d80 100755 (executable)
Binary files a/bin/unres/MD/unres-mult-symetr_ifort_MPICH_GAB.exe and b/bin/unres/MD/unres-mult-symetr_ifort_MPICH_GAB.exe differ

diff --git a/source/unres/src_MD-M/Makefile b/source/unres/src_MD-M/Makefile
deleted file mode 100644 (file)
index 35c2a1f3292e3c562ddb74a675fc20cb0e4cc44b..0000000000000000000000000000000000000000
--- a/source/unres/src_MD-M/Makefile
+++ /dev/null
@@ -1,151 +0,0 @@
-#INSTALL_DIR = /usr/local/mpich-1.2.0
-INSTALL_DIR = /users/software/mpich-1.2.7p1_intel-10.1_em64_ssh
-#
-#FC= /usr/local/opt/intel/compiler60/ia32/bin/ifc
-FC= ifort
-
-OPT =  -O3 -ip -w
-#OPT = -g -CB
-#OPT = -g
-CFLAGS = -DSGI -c
-
-FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include
-FFLAGS1 = -c -w -g -d2 -CA -CB -I$(INSTALL_DIR)/include
-FFLAGS2 = -c -w -g -O0 -I$(INSTALL_DIR)/include
-FFLAGSE = -c -w -O3 -ipo -ipo_obj  -opt_report -I$(INSTALL_DIR)/include
-
-#BIN = ../../../bin/unres/MD-M/unres_Tc_procor_newparm_em64-D-symetr.exe
-#LIBS = -L$(INSTALL_DIR)/lib_pgi -lmpich xdrf/libxdrf.a
-#LIBS = -L$(INSTALL_DIR)/lib_ifort -lmpich xdrf/libxdrf.a
-LIBS = -L$(INSTALL_DIR)/lib -lmpich ../../lib/xdrf_em64/libxdrf.a 
-
-ARCH = LINUX
-PP = /lib/cpp -P
-
-
-all: GAB
-
-.SUFFIXES: .F
-.F.o:
-       ${FC} ${FFLAGS} ${CPPFLAGS} $*.F
-
-object = unres.o arcos.o cartprint.o chainbuild.o convert.o initialize_p.o \
-        matmult.o readrtns_CSA.o parmread.o gen_rand_conf.o printmat.o map.o \
-        pinorm.o randgens.o rescode.o intcor.o timing.o misc.o intlocal.o \
-        cartder.o checkder_p.o econstr_local.o energy_p_new_barrier.o \
-       energy_p_new-sep_barrier.o gradient_p.o minimize_p.o sumsld.o \
-        cored.o rmdd.o geomout.o readpdb.o regularize.o thread.o fitsq.o mcm.o \
-        mc.o bond_move.o refsys.o check_sc_distr.o check_bond.o contact.o djacob.o \
-        eigen.o blas.o add.o entmcm.o \
-        MP.o compare_s1.o \
-        banach.o rmsd.o elecont.o dihed_cons.o \
-        sc_move.o local_move.o \
-        intcartderiv.o lagrangian_lesyng.o\
-       stochfric.o kinetic_lesyng.o MD_A-MTS.o moments.o int_to_cart.o \
-        surfatom.o sort.o muca_md.o rattle.o gauss.o energy_split-sep.o \
-        q_measure.o gnmr1.o test.o ssMD.o permut.o distfit.o checkvar.o
-
-no_option:
-
-GAB: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNRES -DISNAN -DMP -DMPI -DAMD64 \
-       -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC\
-        -DSCCORPDB
-GAB: BIN = ../../../bin/unres/MD-M/unres_ifort_MPICH_GAB.exe
-GAB: ${object} ../../lib/xdrf_em64/libxdrf.a
-       cc -o compinfo compinfo.c
-       ./compinfo | true
-       ${FC} ${FFLAGS} cinfo.f
-       ${FC} ${OPT} ${object} cinfo.o ${LIBS}  -o ${BIN}
-
-E0LL2Y: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNRES -DISNAN -DMP -DMPI -DAMD64 \
-       -DSPLITELE -DLANG0
-E0LL2Y: BIN = ../../../bin/unres/MD-M/unres_ifort_MPICH_E0LL2Y.exe
-E0LL2Y: ${object} ../../lib/xdrf_em64/libxdrf.a
-       cc -o compinfo compinfo.c
-       ./compinfo | true
-       ${FC} ${FFLAGS} cinfo.f
-       ${FC} ${OPT} ${object} cinfo.o ${LIBS}  -o ${BIN}
-
-../../lib/xdrf_em64/libxdrf.a:
-       cd ../../lib/xdrf_em64 && make
-
-4P: CPPFLAGS = -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN \
-       -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC
-4P: BIN = ../../../bin/unres/MD/unres-mult_ifort_single_4P.exe
-4P: ${object} xdrf/libxdrf.a
-       cc -o compinfo compinfo.c
-       ./compinfo | true
-       ${FC} ${FFLAGS} cinfo.f
-       ${FC} ${OPT} ${object} cinfo.o ${LIBS}  -o ${BIN}
-
-E0LL2Y: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN \
-       -DSPLITELE -DLANG0
-E0LL2Y: BIN = ../../../bin/unres/MD/unres-mult_ifort_single_E0LL2Y.exe
-E0LL2Y: ${object} xdrf/libxdrf.a
-       cc -o compinfo compinfo.c
-       ./compinfo | true
-       ${FC} ${FFLAGS} cinfo.f
-       ${FC} ${OPT} ${object} cinfo.o ${LIBS}  -o ${BIN}
-
-xdrf/libxdrf.a:
-       cd xdrf && make
-
-clean:
-       /bin/rm -f *.o && /bin/rm -f compinfo && cd xdrf && make clean
-
-test.o: test.F
-       ${FC} ${FFLAGS} ${CPPFLAGS} test.F
-
-chainbuild.o: chainbuild.F
-       ${FC} ${FFLAGS} ${CPPFLAGS} chainbuild.F
-
-matmult.o: matmult.f
-       ${FC} ${FFLAGS} ${CPPFLAGS} matmult.f
-
-parmread.o : parmread.F
-       ${FC} ${FFLAGS} ${CPPFLAGS} parmread.F
-
-intcor.o : intcor.f
-       ${FC} ${FFLAGS} ${CPPFLAGS} intcor.f
-
-cartder.o : cartder.F
-       ${FC} ${FFLAGS} ${CPPFLAGS} cartder.F
-
-readpdb.o : readpdb.F
-       ${FC} ${FFLAGS2} ${CPPFLAGS} readpdb.F
-
-sumsld.o : sumsld.f
-       ${FC} ${FFLAGS} ${CPPFLAGS} sumsld.f
-        
-cored.o : cored.f
-       ${FC} ${FFLAGS1} ${CPPFLAGS} cored.f
- 
-rmdd.o : rmdd.f
-       ${FC} ${FFLAGS} ${CPPFLAGS} rmdd.f
-
-energy_p_new_barrier.o : energy_p_new_barrier.F
-       ${FC} ${FFLAGSE} ${CPPFLAGS} energy_p_new_barrier.F
-
-gradient_p.o : gradient_p.F
-       ${FC} ${FFLAGSE} ${CPPFLAGS} gradient_p.F
-
-energy_p_new-sep_barrier.o : energy_p_new-sep_barrier.F
-       ${FC} ${FFLAGSE} ${CPPFLAGS} energy_p_new-sep_barrier.F
-
-lagrangian_lesyng.o : lagrangian_lesyng.F
-       ${FC} ${FFLAGSE} ${CPPFLAGS} lagrangian_lesyng.F
-
-MD_A-MTS.o : MD_A-MTS.F
-       ${FC} ${FFLAGSE} ${CPPFLAGS} MD_A-MTS.F
-
-blas.o : blas.f
-       ${FC} ${FFLAGS1} blas.f
-
-add.o : add.f
-       ${FC} ${FFLAGS1} add.f
-
-eigen.o : eigen.f
-       ${FC} ${FFLAGS2} eigen.f
-
-proc_proc.o: proc_proc.c
-       ${CC} ${CFLAGS} proc_proc.c

diff --git a/source/unres/src_MD-M/Makefile b/source/unres/src_MD-M/Makefile
new file mode 120000 (symlink)
index 0000000000000000000000000000000000000000..8453cddeb7b08caf0a886e4bb11a3abc3a14980c
--- /dev/null
+++ b/source/unres/src_MD-M/Makefile
@@ -0,0 +1 @@
+Makefile_MPICH_ifort
\ No newline at end of file

diff --git a/source/unres/src_MD-M/cinfo.f b/source/unres/src_MD-M/cinfo.f
index c7b5141..8757613 100644 (file)
--- a/source/unres/src_MD-M/cinfo.f
+++ b/source/unres/src_MD-M/cinfo.f
@@ -1,15 +1,15 @@
 C DO NOT EDIT THIS FILE - IT HAS BEEN GENERATED BY COMPINFO.C
-C 0 40360 9
+C 0 40360 15
       subroutine cinfo
       include 'COMMON.IOUNITS'
       write(iout,*)'++++ Compile info ++++'
-      write(iout,*)'Version 0.40360 build 9'
-      write(iout,*)'compiled Fri Jan 23 21:00:08 2015'
-      write(iout,*)'compiled by adam@mmka'
+      write(iout,*)'Version 0.40360 build 15'
+      write(iout,*)'compiled Thu Nov 26 21:54:36 2015'
+      write(iout,*)'compiled by adam@piasek4'
       write(iout,*)'OS name:    Linux '
-      write(iout,*)'OS release: 3.2.0-72-generic '
+      write(iout,*)'OS release: 3.2.0-91-generic '
       write(iout,*)'OS version:',
-     & ' #107-Ubuntu SMP Thu Nov 6 14:24:01 UTC 2014 '
+     & ' #129-Ubuntu SMP Wed Sep 9 10:56:06 UTC 2015 '
       write(iout,*)'flags:'
       write(iout,*)'INSTALL_DIR = /users/software/mpich-1.2.7p1_int...'
       write(iout,*)'FC= ifort'

diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F
index c0b4c84..1c6470d 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F
@@ -885,14 +885,16 @@ C Calculate SC interaction energy.
 C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
+c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+c            call flush(iout)
             ind=ind+1
             itypj=itype(j)
             if (itypj.eq.ntyp1) cycle
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
-c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
+c            write (iout,*) "j",j,itypi,itypj,dsc_inv(itypj),dscj_inv,
 c     &       1.0d0/vbld(j+nres)
-c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+c            call flush(iout)
             sig0ij=sigma(itypi,itypj)
             chi1=chi(itypi,itypj)
             chi2=chi(itypj,itypi)
@@ -911,6 +913,8 @@ c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
           yj=mod(yj,boxysize)
           if (yj.lt.0) yj=yj+boxysize
           zj=mod(zj,boxzsize)
+c          write (iout,*) "After box"
+c          call flush(iout)
           if (zj.lt.0) zj=zj+boxzsize
        if ((zj.gt.bordlipbot)
      &.and.(zj.lt.bordliptop)) then
@@ -934,6 +938,8 @@ C lipbufthick is thickenes of lipid buffore
          sslipj=0.0d0
          ssgradlipj=0.0
        endif
+c          write (iout,*) "After lipid"
+c          call flush(iout)
       aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
@@ -981,6 +987,8 @@ C lipbufthick is thickenes of lipid buffore
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
               call sc_angular
+c          write (iout,*) "After sc_angular"
+c          call flush(iout)
               sigsq=1.0D0/sigsq
               sig=sig0ij*dsqrt(sigsq)
               rij_shift=1.0D0/rij-sig+sig0ij
@@ -989,9 +997,9 @@ c              rij_shift=1.2*sig0ij
 C I hate to put IF's in the loops, but here don't have another choice!!!!
               if (rij_shift.le.0.0D0) then
                 evdw=1.0D20
-cd                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-cd     &          restyp(itypi),i,restyp(itypj),j,
-cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
+c                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+c     &          restyp(itypi),i,restyp(itypj),j,
+c     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
                 return
               endif
               sigder=-sig*sigsq
@@ -1005,6 +1013,7 @@ c---------------------------------------------------------------
               eps3der=evdwij*eps2rt
 c              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
+c              call flush(iout)
               evdwij=evdwij*eps2rt*eps3rt
               evdw=evdw+evdwij*sss
               if (lprn) then
@@ -1016,10 +1025,13 @@ c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
      &          eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
      &          evdwij
+                call flush(iout)
               endif
 
-              if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') 
-     &                        'evdw',i,j,evdwij
+              if (energy_dec) then
+                write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij
+                call flush(iout)
+              endif
 
 C Calculate gradient components.
               e1=e1*eps1*eps2rt**2*eps3rt**2
@@ -1034,8 +1046,12 @@ C Calculate the radial part of the gradient
               gg(3)=zj*fac
               gg_lipi(3)=ssgradlipi*evdwij
               gg_lipj(3)=ssgradlipj*evdwij
+c              write (iout,*) "Calling sc_grad_scale"
+c              call flush(iout)
 C Calculate angular part of the gradient.
               call sc_grad_scale(sss)
+c              write (iout,*) "After sc_grad_scale"
+c              call flush(iout)
             endif
           enddo      ! j
         enddo        ! iint
@@ -1298,6 +1314,8 @@ C----------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       double precision dcosom1(3),dcosom2(3)
       double precision scalfac
+c      write (iout,*) "sc_grad_scale",i,j,k,l
+c      call flush(iout)
       eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
       eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
       eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
@@ -1335,8 +1353,8 @@ C
 C Calculate the components of the gradient in DC and X
 C
       do l=1,3
-        gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k)
-        gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k)
+        gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
+        gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
       enddo
       return
       end

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index 00862f2..747b132 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -130,6 +130,8 @@ cmc
 c      if (dyn_ss) call dyn_set_nss
 
 c      print *,"Processor",myrank," computed USCSC"
+c      write (iout,*) "SCSC computed OK"
+c      call flush_(iout)
 #ifdef TIMING
       time01=MPI_Wtime() 
 #endif
@@ -171,6 +173,8 @@ c      print *,"Processor",myrank," left VEC_AND_DERIV"
 c        call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
 c     &   eello_turn4)
       endif
+c      write (iout,*) "eelec computed OK"
+c      call flush_(iout)
 c      print *,"Processor",myrank," computed UELEC"
 C
 C Calculate excluded-volume interaction energy between peptide groups
@@ -187,10 +191,14 @@ C
 c        write (iout,*) "Soft-sphere SCP potential"
         call escp_soft_sphere(evdw2,evdw2_14)
       endif
+c      write (iout,*) "escp computed OK"
+c      call flush_(iout)
 c
 c Calculate the bond-stretching energy
 c
       call ebond(estr)
+c      write (iout,*) "ebond computed OK"
+c      call flush_(iout)
 C 
 C Calculate the disulfide-bridge and other energy and the contributions
 C from other distance constraints.
@@ -206,12 +214,16 @@ C
         ebe=0
         ethetacnstr=0
       endif
+c      write (iout,*) "ebend computed OK"
+c      call flush_(iout)
 c      print *,"Processor",myrank," computed UB"
 C
 C Calculate the SC local energy.
 C
 C      print *,"TU DOCHODZE?"
       call esc(escloc)
+c      write (iout,*) "esc computed OK"
+c      call flush_(iout)
 c      print *,"Processor",myrank," computed USC"
 C
 C Calculate the virtual-bond torsional energy.
@@ -223,6 +235,8 @@ cd    print *,'nterm=',nterm
        etors=0
        edihcnstr=0
       endif
+c      write (iout,*) "etor computed OK"
+c      call flush_(iout)
 c      print *,"Processor",myrank," computed Utor"
 C
 C 6/23/01 Calculate double-torsional energy
@@ -232,6 +246,8 @@ C
       else
        etors_d=0
       endif
+c      write (iout,*) "etor_d computed OK"
+c      call flush_(iout)
 c      print *,"Processor",myrank," computed Utord"
 C
 C 21/5/07 Calculate local sicdechain correlation energy
@@ -241,6 +257,8 @@ C
       else
         esccor=0.0d0
       endif
+c      write (iout,*) "eback_sc_corr computed OK"
+c      call flush_(iout)
 C      print *,"PRZED MULIt"
 c      print *,"Processor",myrank," computed Usccorr"
 C 
@@ -250,9 +268,13 @@ C
       n_corr1=0
       if ((wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 
      &    .or. wturn6.gt.0.0d0) .and. ipot.lt.6) then
+c         write (iout,*) "Calling multibody_eello"
+c         call flush_(iout)
          call multibody_eello(ecorr,ecorr5,ecorr6,eturn6,n_corr,n_corr1)
-cd         write(2,*)'multibody_eello n_corr=',n_corr,' n_corr1=',n_corr1,
-cd     &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6
+c         write(iout,*)
+c     & 'multibody_eello n_corr=',n_corr,' n_corr1=',n_corr1,
+c     & " ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6
+c         call flush_(iout)
       else
          ecorr=0.0d0
          ecorr5=0.0d0
@@ -260,9 +282,14 @@ cd     &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6
          eturn6=0.0d0
       endif
       if ((wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) .and. ipot.lt.6) then
+c         write (iout,*) "Calling multibody_gb_ecorr"
+c         call flush_(iout)
          call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
-cd         write (iout,*) "multibody_hb ecorr",ecorr
+c         write (iout,*) "Exited multibody_hb ecorr",ecorr
+c         call flush_(iout)
       endif
+c      write (iout,*) "multibody computed OK"
+c      call flush_(iout)
 c      print *,"Processor",myrank," computed Ucorr"
 C 
 C If performing constraint dynamics, call the constraint energy
@@ -281,12 +308,15 @@ C      print *,"przed lipidami"
       if (wliptran.gt.0) then
         call Eliptransfer(eliptran)
       endif
-C      print *,"za lipidami"
+c      write (iout,*) "lipid energy computed OK"
+c      call flush_(iout)
       if (AFMlog.gt.0) then
         call AFMforce(Eafmforce)
       else if (selfguide.gt.0) then
         call AFMvel(Eafmforce)
       endif
+c      write (iout,*) "AFMforce computed OK"
+c      call flush_(iout)
 #ifdef TIMING
       time_enecalc=time_enecalc+MPI_Wtime()-time00
 #endif
@@ -335,7 +365,11 @@ c    Here are the energies showed per procesor if the are more processors
 c    per molecule then we sum it up in sum_energy subroutine 
 c      print *," Processor",myrank," calls SUM_ENERGY"
       call sum_energy(energia,.true.)
+c      write (iout,*) "sum energy OK"
+c      call flush_(iout)
       if (dyn_ss) call dyn_set_nss
+c      write (iout,*) "Exiting energy"
+c      call flush_(iout)
 c      print *," Processor",myrank," left SUM_ENERGY"
 #ifdef TIMING
       time_sumene=time_sumene+MPI_Wtime()-time00
@@ -5638,6 +5672,8 @@ c      time12=1.0d0
       etheta=0.0D0
 c     write (*,'(a,i2)') 'EBEND ICG=',icg
       do i=ithet_start,ithet_end
+c        write (iout,*) "ebend: i=",i
+c        call flush_(iout)
         if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
      &  .or.itype(i).eq.ntyp1) cycle
 C Zero the energy function and its derivative at 0 or pi.
@@ -5739,9 +5775,12 @@ C Derivatives of the "mean" values in gamma1 and gamma2.
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
         gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg)
       enddo
+c      write (iout,*) "Exit loop" 
+c      call flush_(iout)
       ethetacnstr=0.0d0
-C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
-      do i=ithetaconstr_start,ithetaconstr_end
+c      write (iout,*) ithetaconstr_start,ithetaconstr_end,"TU"
+c      call flush_(iout)
+      do i=max0(ithetaconstr_start,1),ithetaconstr_end
         itheta=itheta_constr(i)
         thetiii=theta(itheta)
         difi=pinorm(thetiii-theta_constr0(i))
@@ -5766,6 +5805,8 @@ C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
      &    gloc(itheta+nphi-2,icg)
         endif
       enddo
+c      write (iout,*) "Exit ebend"
+c      call flush_(iout)
 
 C Ufff.... We've done all this!!! 
       return
@@ -6103,7 +6144,7 @@ c        lprn1=.false.
 C now constrains
       ethetacnstr=0.0d0
 C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
-      do i=ithetaconstr_start,ithetaconstr_end
+      do i=max0(ithetaconstr_start,1),ithetaconstr_end
         itheta=itheta_constr(i)
         thetiii=theta(itheta)
         difi=pinorm(thetiii-theta_constr0(i))
@@ -6986,7 +7027,7 @@ c----------------------------------------------------------------------------
       logical lprn
 C Set lprn=.true. for debugging
       lprn=.false.
-c     lprn=.true.
+c      lprn=.true.
       etors=0.0D0
       do i=iphi_start,iphi_end
 C ANY TWO ARE DUMMY ATOMS in row CYCLE
@@ -7379,7 +7420,7 @@ C Set lprn=.true. for debugging
       if (lprn) then
         write (iout,'(a)') 'Contact function values before RECEIVE:'
         do i=nnt,nct-2
-          write (iout,'(2i3,50(1x,i2,f5.2))') 
+          write (iout,'(2i3,50(1x,i3,f5.2))') 
      &    i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i),
      &    j=1,num_cont_hb(i))
         enddo
@@ -7618,6 +7659,7 @@ C Calculate the local-electrostatic correlation terms
             jp1=iabs(j1)
 c            write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
 c     &         ' jj=',jj,' kk=',kk
+c            call flush(iout)
             if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0 
      &          .or. j.lt.0 .and. j1.gt.0) .and.
      &         (jp1.eq.jp+1 .or. jp1.eq.jp-1)) then
@@ -7635,8 +7677,9 @@ c             ecorr=ecorr+ehbcorr(i,j,i+1,j,jj,kk,0.60D0,-0.40D0)
           enddo ! kk
           do kk=1,num_conti
             j1=jcont_hb(kk,i)
-c           write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
-c    &         ' jj=',jj,' kk=',kk
+c            write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
+c     &         ' jj=',jj,' kk=',kk
+c            call flush(iout)
             if (j1.eq.j+1) then
 C Contacts I-J and (I+1)-J occur simultaneously. 
 C The system loses extra energy.
@@ -7730,6 +7773,7 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding
       include 'COMMON.CONTACTS'
       include 'COMMON.CHAIN'
       include 'COMMON.CONTROL'
+      include 'COMMON.TORSION'
       double precision gx(3),gx1(3)
       integer num_cont_hb_old(maxres)
       logical lprn,ldone
@@ -7738,6 +7782,8 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding
 C Set lprn=.true. for debugging
       lprn=.false.
       eturn6=0.0d0
+c      write (iout,*) "MULTIBODY_EELLO"
+c      call flush(iout)
 #ifdef MPI
       do i=1,nres
         num_cont_hb_old(i)=num_cont_hb(i)
@@ -7748,12 +7794,12 @@ C Set lprn=.true. for debugging
       if (lprn) then
         write (iout,'(a)') 'Contact function values before RECEIVE:'
         do i=nnt,nct-2
-          write (iout,'(2i3,50(1x,i2,f5.2))') 
+          write (iout,'(2i3,50(1x,i3,f5.2))') 
      &    i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i),
      &    j=1,num_cont_hb(i))
         enddo
+        call flush(iout)
       endif
-      call flush(iout)
       do i=1,ntask_cont_from
         ncont_recv(i)=0
       enddo
@@ -7955,10 +8001,15 @@ c          call flush(iout)
       if (lprn) then
         write (iout,'(a)') 'Contact function values:'
         do i=nnt,nct-2
-          write (iout,'(2i3,50(1x,i2,5f6.3))') 
+          write (iout,'(2i3,50(1x,i3,5f6.3))') 
      &    i,num_cont_hb(i),(jcont_hb(j,i),d_cont(j,i),
      &    ((a_chuj(ll,kk,j,i),ll=1,2),kk=1,2),j=1,num_cont_hb(i))
         enddo
+        write (iout,*) "itortyp"
+        do i=1,nres
+          write (iout,*) i,itype(i),itortyp(itype(i))
+        enddo
+        call flush(iout)
       endif
       ecorr=0.0D0
       ecorr5=0.0d0
@@ -7999,8 +8050,11 @@ c        write (iout,*) "corr loop i",i
           do kk=1,num_conti1
             j1=jcont_hb(kk,i1)
             jp1=iabs(j1)
-c            write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
-c     &         ' jj=',jj,' kk=',kk
+            if (lprn) then
+              write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
+     &         ' jj=',jj,' kk=',kk
+              call flush(iout)
+            endif
 c            if (j1.eq.j+1 .or. j1.eq.j-1) then
             if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0 
      &          .or. j.lt.0 .and. j1.gt.0) .and.
@@ -8040,6 +8094,12 @@ c                do iii=1,nres
 c                  write (iout,'(i5,3f10.5)') 
 c     &             iii,(gradcorr5(jjj,iii),jjj=1,3)
 c                enddo
+c                write (iout,*) "ecorr4"
+c                call flush(iout)
+c                write (iout,*) "eello5:",i,jp,i+1,jp1,jj,kk,
+c     &        itype(jp),itype(i+1),itype(jp1),
+c     &        itortyp(itype(jp)),itortyp(itype(i+1)),itortyp(itype(jp1))
+c                call flush(iout)
                 if (wcorr5.gt.0.0d0)
      &            ecorr5=ecorr5+eello5(i,jp,i+1,jp1,jj,kk)
 c                write (iout,*) "gradcorr5 after eello5"
@@ -8052,12 +8112,16 @@ c                enddo
      2                'ecorr5',i,j,i+1,j1,eello5(i,jp,i+1,jp1,jj,kk)
 cd                write(2,*)'wcorr6',wcorr6,' wturn6',wturn6
 cd                write(2,*)'ijkl',i,jp,i+1,jp1 
+c                write (iout,*) "ecorr5"
+c                call flush(iout)
                 if (wcorr6.gt.0.0d0 .and. (jp.ne.i+4 .or. jp1.ne.i+3
      &               .or. wturn6.eq.0.0d0))then
 cd                  write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1
                   ecorr6=ecorr6+eello6(i,jp,i+1,jp1,jj,kk)
                   if (energy_dec) write (iout,'(a6,4i5,0pf7.3)')
      1                'ecorr6',i,j,i+1,j1,eello6(i,jp,i+1,jp1,jj,kk)
+c                write (iout,*) "ecorr6"
+c                call flush(iout)
 cd                write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5,
 cd     &            'ecorr6=',ecorr6
 cd                write (iout,'(4e15.5)') sred_geom,
@@ -8071,10 +8135,13 @@ cd                  write (iout,*) '******eturn6: i,j,i+1,j1',i,jip,i+1,jp1
                   if (energy_dec) write (iout,'(a6,4i5,0pf7.3)')
      1                 'eturn6',i,j,i+1,j1,eello_turn6(i,jj,kk)
 cd                  write (2,*) 'multibody_eello:eturn6',eturn6
+c                write (iout,*) "ecorr4"
+c                call flush(iout)
                 endif
               ENDIF
 1111          continue
             endif
+                  if (energy_dec) call flush(iout)
           enddo ! kk
         enddo ! jj
       enddo ! i
@@ -8837,9 +8904,9 @@ cd      endif
 cd      write (iout,*)
 cd     &   'EELLO5: Contacts have occurred for peptide groups',i,j,
 cd     &   ' and',k,l
-      itk=itortyp(itype(k))
-      itl=itortyp(itype(l))
-      itj=itortyp(itype(j))
+c      itk=itortyp(itype(k))
+c      itl=itortyp(itype(l))
+c      itj=itortyp(itype(j))
       eello5_1=0.0d0
       eello5_2=0.0d0
       eello5_3=0.0d0
@@ -8908,7 +8975,7 @@ C Cartesian gradient
 c      goto 1112
 c1111  continue
 C Contribution from graph II 
-      call transpose2(EE(1,1,itk),auxmat(1,1))
+      call transpose2(EE(1,1,k),auxmat(1,1))
       call matmat2(auxmat(1,1),AEA(1,1,1),pizda(1,1))
       vv(1)=pizda(1,1)+pizda(2,2)
       vv(2)=pizda(2,1)-pizda(1,2)
@@ -8989,7 +9056,7 @@ C Cartesian gradient
 cd        goto 1112
 C Contribution from graph IV
 cd1110    continue
-        call transpose2(EE(1,1,itl),auxmat(1,1))
+        call transpose2(EE(1,1,l),auxmat(1,1))
         call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
         vv(1)=pizda(1,1)+pizda(2,2)
         vv(2)=pizda(2,1)-pizda(1,2)
@@ -9062,7 +9129,7 @@ C Cartesian gradient
 cd        goto 1112
 C Contribution from graph IV
 1110    continue
-        call transpose2(EE(1,1,itj),auxmat(1,1))
+        call transpose2(EE(1,1,j),auxmat(1,1))
         call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
         vv(1)=pizda(1,1)+pizda(2,2)
         vv(2)=pizda(2,1)-pizda(1,2)
@@ -9359,7 +9426,7 @@ C       o     o       o     o                                                  C
 C       i             i                                                        C
 C                                                                              C
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-      itk=itortyp(itype(k))
+c      itk=itortyp(itype(k))
       s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i))
       s2=-scalar2(AEAb2(1,1,imat),Ug2Db1t(1,k))
       s3= scalar2(AEAb2(1,1,imat),CUgb2(1,k))
@@ -9649,19 +9716,19 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective 
 C           energy moment and not to the cluster cumulant.
-      iti=itortyp(itype(i))
-      if (j.lt.nres-1) then
-        itj1=itortyp(itype(j+1))
-      else
-        itj1=ntortyp
-      endif
-      itk=itortyp(itype(k))
-      itk1=itortyp(itype(k+1))
-      if (l.lt.nres-1) then
-        itl1=itortyp(itype(l+1))
-      else
-        itl1=ntortyp
-      endif
+c      iti=itortyp(itype(i))
+c      if (j.lt.nres-1) then
+c        itj1=itortyp(itype(j+1))
+c      else
+c        itj1=ntortyp
+c      endif
+c      itk=itortyp(itype(k))
+c      itk1=itortyp(itype(k+1))
+c      if (l.lt.nres-1) then
+c        itl1=itortyp(itype(l+1))
+c      else
+c        itl1=ntortyp
+c      endif
 #ifdef MOMENT
       s1=dip(4,jj,i)*dip(4,kk,k)
 #endif
@@ -9669,7 +9736,7 @@ C           energy moment and not to the cluster cumulant.
       s2=0.5d0*scalar2(b1(1,k),auxvec(1))
       call matvec2(AECA(1,1,2),b1(1,l+1),auxvec(1))
       s3=0.5d0*scalar2(b1(1,j+1),auxvec(1))
-      call transpose2(EE(1,1,itk),auxmat(1,1))
+      call transpose2(EE(1,1,k),auxmat(1,1))
       call matmat2(auxmat(1,1),AECA(1,1,1),pizda(1,1))
       vv(1)=pizda(1,1)+pizda(2,2)
       vv(2)=pizda(2,1)-pizda(1,2)
@@ -9767,25 +9834,25 @@ C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective 
 C           energy moment and not to the cluster cumulant.
 cd      write (2,*) 'eello_graph4: wturn6',wturn6
-      iti=itortyp(itype(i))
-      itj=itortyp(itype(j))
-      if (j.lt.nres-1) then
-        itj1=itortyp(itype(j+1))
-      else
-        itj1=ntortyp
-      endif
-      itk=itortyp(itype(k))
-      if (k.lt.nres-1) then
-        itk1=itortyp(itype(k+1))
-      else
-        itk1=ntortyp
-      endif
-      itl=itortyp(itype(l))
-      if (l.lt.nres-1) then
-        itl1=itortyp(itype(l+1))
-      else
-        itl1=ntortyp
-      endif
+c      iti=itortyp(itype(i))
+c      itj=itortyp(itype(j))
+c      if (j.lt.nres-1) then
+c        itj1=itortyp(itype(j+1))
+c      else
+c        itj1=ntortyp
+c      endif
+c      itk=itortyp(itype(k))
+c      if (k.lt.nres-1) then
+c        itk1=itortyp(itype(k+1))
+c      else
+c        itk1=ntortyp
+c      endif
+c      itl=itortyp(itype(l))
+c      if (l.lt.nres-1) then
+c        itl1=itortyp(itype(l+1))
+c      else
+c        itl1=ntortyp
+c      endif
 cd      write (2,*) 'eello6_graph4:','i',i,' j',j,' k',k,' l',l
 cd      write (2,*) 'iti',iti,' itj',itj,' itj1',itj1,' itk',itk,
 cd     & ' itl',itl,' itl1',itl1
@@ -10798,7 +10865,7 @@ C  gradient po costhet
      &*VofOverlap
 C grad_shield_side is Cbeta sidechain gradient
       grad_shield_side(j,ishield_list(i),i)=
-     &        (sh_frac_dist_grad(j)*-2.0d0
+     &        (sh_frac_dist_grad(j)*(-2.0d0)
      &       +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet)
      &       +scale_fac_dist*(cosphi_grad_long(j))
      &        *2.0d0/(1.0-cosphi))

diff --git a/source/unres/src_MD-M/energy_split-sep.F b/source/unres/src_MD-M/energy_split-sep.F
index 24ab8dd..b9fc9e0 100644 (file)
--- a/source/unres/src_MD-M/energy_split-sep.F
+++ b/source/unres/src_MD-M/energy_split-sep.F
@@ -36,7 +36,7 @@ c        if (fg_rank.eq.0) call int_from_cart1(.false.)
 #ifdef MPI      
 c      write(iout,*) "ETOTAL_LONG Processor",fg_rank,
 c     & " absolute rank",myrank," nfgtasks",nfgtasks
-      call flush(iout)
+c      call flush(iout)
       if (nfgtasks.gt.1) then
         time00=MPI_Wtime()
 C FG slaves call the following matching MPI_Bcast in ERGASTULUM
@@ -369,6 +369,8 @@ c      call int_from_cart1(.false.)
 C
 C Compute the side-chain and electrostatic interaction energy
 C
+c      write (iout,*) "Computing SC"
+c      call flush(iout)
       goto (101,102,103,104,105,106) ipot
 C Lennard-Jones potential.
   101 call elj_short(evdw)
@@ -396,43 +398,66 @@ C
 c
 c Calculate the short-range part of Evdwpp
 c
+c      write (iout,*) "SC_short computed OK"
+c      call flush(iout)
       call evdwpp_short(evdw1)
+c      write (iout,*) "VDWPP_short computed OK"
+c      call flush(iout)
 c
 c Calculate the short-range part of ESCp
 c
       if (ipot.lt.6) then
         call escp_short(evdw2,evdw2_14)
+c      write (iout,*) "SCP_short computed OK"
+c      call flush(iout)
       endif
 c
 c Calculate the bond-stretching energy
 c
       call ebond(estr)
+c      write (iout,*) "ebond_short computed OK"
+c      call flush(iout)
 C 
 C Calculate the disulfide-bridge and other energy and the contributions
 C from other distance constraints.
       call edis(ehpb)
+c      write (iout,*) "ehpb_short computed OK"
+c      call flush(iout)
 C
 C Calculate the virtual-bond-angle energy.
 C
-      call ebend(ebe)
+      ethetacnstr=0.0d0
+      call ebend(ebe,ethetacnstr)
+c      write (iout,*) "ebe_short computed OK"
+c      call flush(iout)
 C
 C Calculate the SC local energy.
 C
       call vec_and_deriv
+c      write (iout,*) "vec_and_dervi_short computed OK"
+c      call flush(iout)
       call esc(escloc)
+c      write (iout,*) "SCLOCAL_short computed OK"
+c      call flush(iout)
 C
 C Calculate the virtual-bond torsional energy.
 C
       call etor(etors,edihcnstr)
+c      write (iout,*) "TOR_short computed OK"
+c      call flush(iout)
 C
 C 6/23/01 Calculate double-torsional energy
 C
       call etor_d(etors_d)
+c      write (iout,*) "DTOR_short computed OK"
+c      call flush(iout)
 C
 C 21/5/07 Calculate local sicdechain correlation energy
 C
       if (wsccor.gt.0.0d0) then
         call eback_sc_corr(esccor)
+c      write (iout,*) "SCCORR_short computed OK"
+c      call flush(iout)
       else
         esccor=0.0d0
       endif
@@ -464,9 +489,9 @@ C
       energia(19)=edihcnstr
       energia(21)=esccor
 c      write (iout,*) "ETOTAL_SHORT before SUM_ENERGY"
-      call flush(iout)
+c      call flush(iout)
       call sum_energy(energia,.true.)
 c      write (iout,*) "Exit ETOTAL_SHORT"
-      call flush(iout)
+c      call flush(iout)
       return
       end

diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F
index ae4d710..23fe7df 100644 (file)
--- a/source/unres/src_MD-M/parmread.F
+++ b/source/unres/src_MD-M/parmread.F
@@ -58,7 +58,7 @@ C Assign virtual-bond length
       vblinv2=vblinv*vblinv
 c
 c Read the virtual-bond parameters, masses, and moments of inertia
-c and Stokes' radii of the peptide group and side chains
+c and Stokes radii of the peptide group and side chains
 c
 #ifdef CRYST_BOND
       read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok
@@ -99,13 +99,30 @@ c
      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
           enddo
         enddo
+#ifdef AIX
+        call flush_(iout)
+#else
+        call flush(iout)
+#endif
       endif
 C reading lipid parameters
-       read(iliptranpar,*) pepliptran
+       read(iliptranpar,*,end=120,err=120) pepliptran
        do i=1,ntyp
-       read(iliptranpar,*) liptranene(i)
+       read(iliptranpar,*,end=120,err=120) liptranene(i)
        enddo
        close(iliptranpar)
+       if (lprint) then
+         write (iout,*) "Lipid transfer parameters"
+         write (iout,'(a5,f10.5)') "pept",pepliptran
+         do i=1,ntyp
+           write (iout,'(a5,f10.5)') restyp(i),liptranene(i)
+         enddo
+#ifdef AIX
+        call flush_(iout)
+#else
+        call flush(iout)
+#endif
+       endif
 #ifdef CRYST_THETA
 C
 C Read the parameters of the probability distribution/energy expression 
@@ -182,6 +199,7 @@ C
           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
      &       sig0(i),(gthet(j,i),j=1,3)
         enddo
+        call flush(iout)
        else
        write (iout,'(a)') 
      &    'Parameters of the virtual-bond valence angles:'
@@ -210,6 +228,11 @@ C
           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
        enddo
+#ifdef AIX
+        call flush_(iout)
+#else
+        call flush(iout)
+#endif
       endif
       endif
 #else 
@@ -395,9 +418,13 @@ C
           enddo
         enddo
       enddo
-      call flush(iout)
+#ifdef AIX
+        call flush_(iout)
+#else
+        call flush(iout)
+#endif
       endif
-      write (2,*) "Start reading THETA_PDB",ithep_pdb
+      write (iout,*) "Start reading THETA_PDB",ithep_pdb
       do i=1,ntyp
 c      write (2,*) 'i=',i
         read (ithep_pdb,*,err=111,end=111)
@@ -442,7 +469,7 @@ c      write (2,*) 'i=',i
          gthet(j,i)=gthet(j,-i)
        enddo
       enddo
-      write (2,*) "End reading THETA_PDB"
+      write (iout,*) "End reading THETA_PDB"
       close (ithep_pdb)
 #endif
       close(ithep)
@@ -530,6 +557,11 @@ C BSC is amplitude of Gaussian
             endif
          endif
         enddo
+#ifdef AIX
+        call flush_(iout)
+#else
+        call flush(iout)
+#endif
       endif
 #else
 C 
@@ -550,7 +582,7 @@ C
 C Read the parameters of the probability distribution/energy expression
 C of the side chains.
 C
-      write (2,*) "Start reading ROTAM_PDB"
+      write (iout,*) "Start reading ROTAM_PDB"
       do i=1,ntyp
         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
         if (i.eq.10) then
@@ -589,7 +621,12 @@ C
         endif
       enddo
       close (irotam_pdb)
-      write (2,*) "End reading ROTAM_PDB"
+c      write (iout,*) "End reading ROTAM_PDB"
+c#ifdef AIX
+c      call flush_(iout)
+c#else
+c      call flush(iout)
+c#endif
 #endif
       close(irotam)
 
@@ -617,6 +654,11 @@ C
            write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
           enddo
         enddo
+#ifdef AIX
+        call flush_(iout)
+#else
+        call flush(iout)
+#endif
       endif
 #else
 C
@@ -659,8 +701,16 @@ c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
       enddo
       enddo
       close (itorp)
+c      write (iout,*) "End reading torsional parameters"
+c#ifdef AIX
+c      call flush_(iout)
+c#else
+c      call flush(iout)
+c#endif
       if (lprint) then
         write (iout,'(/a/)') 'Torsional constants:'
+        do iblock=1,2
+        write (iout,*) "IBLOCK",iblock
         do i=1,ntortyp
           do j=1,ntortyp
             write (iout,*) 'ityp',i,' jtyp',j
@@ -676,6 +726,12 @@ c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
             enddo
           enddo
         enddo
+        enddo
+#ifdef AIX
+        call flush_(iout)
+#else
+        call flush(iout)
+#endif
       endif
 
 C
@@ -773,6 +829,11 @@ C Martix of D parameters for two dimesional fourier series
         enddo
       enddo
       enddo
+#ifdef AIX
+        call flush_(iout)
+#else
+        call flush(iout)
+#endif
       endif
 #endif
 C Read of Side-chain backbone correlation parameters
@@ -780,6 +841,12 @@ C Modified 11 May 2012 by Adasko
 CCC
 C
       read (isccor,*,end=119,err=119) nsccortyp
+c      write (iout,*) "Reading sccor parameters",nsccortyp
+c#ifdef AIX
+c      call flush_(iout)
+c#else
+c      call flush(iout)
+c#endif
 #ifdef SCCORPDB
       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
       do i=-ntyp,-1
@@ -882,15 +949,23 @@ cc maxinter is maximum interaction sites
             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
      &(1+vlor3sccor(k,i,j)**2)
           enddo
-          v0sccor(i,j,iblock)=v0ijsccor
+          v0sccor(l,i,j)=v0ijsccor
         enddo
       enddo
       enddo
       close (isccor)
 
 #endif      
+c      write (iout,*) "sccor parameters read"
+c#ifdef AIX
+c      call flush_(iout)
+c#else
+c      call flush(iout)
+c#endif
       if (lprint) then
-        write (iout,'(/a/)') 'Torsional constants:'
+        write (iout,'(/a/)') 'SC-torsional constants:'
+        do l=1,maxinter
+        write (iout,*) "Torsional type",l
         do i=1,nsccortyp
           do j=1,nsccortyp
             write (iout,*) 'ityp',i,' jtyp',j
@@ -905,6 +980,12 @@ cc maxinter is maximum interaction sites
             enddo
           enddo
         enddo
+        enddo
+#ifdef AIX
+        call flush_(iout)
+#else
+        call flush(iout)
+#endif
       endif
 
 C
@@ -1047,6 +1128,11 @@ c      lprint=.true.
           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
         enddo
       enddo
+#ifdef AIX
+        call flush_(iout)
+#else
+        call flush(iout)
+#endif
       endif
 c      lprint=.false.
 
@@ -1077,6 +1163,11 @@ c        lprint=.true.
 c        lprint=.false.
         enddo
       enddo
+#ifdef AIX
+        if (lprint) call flush_(iout)
+#else
+        if (lprint) call flush(iout)
+#endif
 C
 C Read side-chain interaction parameters.
 C
@@ -1251,6 +1342,11 @@ c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
          endif
         enddo
       enddo
+#ifdef AIX
+        if (lprint) call flush_(iout)
+#else
+        if (lprint) call flush(iout)
+#endif
 #ifdef OLDSCP
 C
 C Define the SC-p interaction constants (hard-coded; old style)
@@ -1297,6 +1393,11 @@ c      lprint=.true.
           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
         enddo
+#ifdef AIX
+        call flush_(iout)
+#else
+        call flush(iout)
+#endif
       endif
 c      lprint=.false.
 #endif
@@ -1357,25 +1458,82 @@ C      buff_shield=1.0d0
 C      endif
       return
   111 write (iout,*) "Error reading bending energy parameters."
+#ifdef AIX
+      call flush_(iout)
+#else
+      call flush(iout)
+#endif
       goto 999
   112 write (iout,*) "Error reading rotamer energy parameters."
+#ifdef AIX
+      call flush_(iout)
+#else
+      call flush(iout)
+#endif
       goto 999
   113 write (iout,*) "Error reading torsional energy parameters."
+#ifdef AIX
+      call flush_(iout)
+#else
+      call flush(iout)
+#endif
       goto 999
   114 write (iout,*) "Error reading double torsional energy parameters."
+#ifdef AIX
+      call flush_(iout)
+#else
+      call flush(iout)
+#endif
       goto 999
   115 write (iout,*) 
      &  "Error reading cumulant (multibody energy) parameters."
+#ifdef AIX
+      call flush_(iout)
+#else
+      call flush(iout)
+#endif
       goto 999
   116 write (iout,*) "Error reading electrostatic energy parameters."
+#ifdef AIX
+      call flush_(iout)
+#else
+      call flush(iout)
+#endif
       goto 999
  1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
+#ifdef AIX
+      call flush_(iout)
+#else
+      call flush(iout)
+#endif
       goto 999
   117 write (iout,*) "Error reading side chain interaction parameters."
+#ifdef AIX
+      call flush_(iout)
+#else
+      call flush(iout)
+#endif
       goto 999
   118 write (iout,*) "Error reading SCp interaction parameters."
+#ifdef AIX
+      call flush_(iout)
+#else
+      call flush(iout)
+#endif
       goto 999
   119 write (iout,*) "Error reading SCCOR parameters"
+#ifdef AIX
+      call flush_(iout)
+#else
+      call flush(iout)
+#endif
+      goto 999
+  120 write (iout,*) "Error reading lipid parameters"
+#ifdef AIX
+      call flush_(iout)
+#else
+      call flush(iout)
+#endif
   999 continue
 #ifdef MPI
       call MPI_Finalize(Ierror)