Merge branch 'devel' into AFM
authorAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Thu, 6 Aug 2015 08:39:41 +0000 (10:39 +0200)
committerAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Thu, 6 Aug 2015 08:39:41 +0000 (10:39 +0200)
Conflicts:
.gitignore
PARAM/pot_theta_G631_DIL.parm
bin/unres/MD/unres-mult-symetr_ifort_MPICH_E0LL2Y.exe
bin/unres/MD/unres_Tc_procor_oldparm_em64-D-symetr.exe
bin/unres/MD/unres_ifort_MPICH_GAB.exe
bin/xdrf2ang
bin/xdrf2pdb
bin/xdrf2pdb-m
source/cluster/wham/src-M/energy_p_new.F
source/maxlik/src_CSA/CMakeLists.txt
source/unres/src_CSA_DiL/CMakeLists.txt
source/unres/src_CSA_DiL/Makefile
source/unres/src_CSA_DiL/csa.F
source/unres/src_MD-M/COMMON.CONTROL
source/unres/src_MD-M/DIMENSIONS
source/unres/src_MD-M/Makefile
source/unres/src_MD-M/cinfo.f
source/unres/src_MD-M/energy_p_new_barrier.F
source/unres/src_MD-M/parmread.F
source/unres/src_MD-M/refsys.f
source/unres/src_MD/COMMON.TORSION
source/wham/src-M/DIMENSIONS
source/wham/src-M/energy_p_new.F
source/wham/src-M/readrtns.F
source/xdrfpdb/src/Makefile

29 files changed:
1  2 
.gitignore
bin/unres/MD/unres-mult-symetr_ifort_MPICH_E0LL2Y.exe
bin/unres/MD/unres_ifort_MPICH_GAB.exe
source/cluster/wham/src-M/energy_p_new.F
source/cluster/wham/src-M/parmread.F
source/cluster/wham/src-M/readrtns.F
source/unres/src_CSA/CMakeLists.txt
source/unres/src_MD-M/COMMON.CONTROL
source/unres/src_MD-M/DIMENSIONS
source/unres/src_MD-M/Makefile
source/unres/src_MD-M/checkder_p.F
source/unres/src_MD-M/energy_p_new_barrier.F
source/unres/src_MD-M/geomout.F
source/unres/src_MD-M/initialize_p.F
source/unres/src_MD-M/intcartderiv.F
source/unres/src_MD-M/parmread.F
source/unres/src_MD-M/readpdb.F
source/unres/src_MD-M/readrtns_CSA.F
source/unres/src_MD-M/refsys.f
source/unres/src_MD-M/ssMD.F
source/unres/src_MD-M/unres.F
source/unres/src_MD/initialize_p.F
source/unres/src_MD/parmread.F
source/unres/src_MD/ssMD.F
source/wham/src-M/enecalc1.F
source/wham/src-M/energy_p_new.F
source/wham/src-M/initialize_p.F
source/wham/src-M/parmread.F
source/wham/src-M/readrtns.F

diff --cc .gitignore
@@@ -17,10 -16,6 +16,7 @@@ build_prere
  period_build/
  period_build2/
  build_*/
- =======
- build_*/
 +period_*/
- >>>>>>> master
  
  # latex files in documentation 
  doc/*/*.aux
diff --cc bin/unres/MD/unres-mult-symetr_ifort_MPICH_E0LL2Y.exe
index a440770,3128072..0000000
deleted file mode 100755,100755
Binary files differ
diff --cc bin/unres/MD/unres_ifort_MPICH_GAB.exe
index 2f98745,c31777f..0000000
deleted file mode 100755,100755
Binary files differ
        common /srutu/icall
        integer icant
        external icant
 -      logical energy_dec /.true./
 +      integer xshift,yshift,zshift
++      logical energy_dec /.false./
  c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=0.0D0
        evdw_t=0.0d0
@@@ -3326,15 -3280,11 +3461,15 @@@ C Zero the energy function and its deri
            ichir21=isign(1,itype(i))
            ichir22=isign(1,itype(i))
           endif
 -        if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
 +         if (i.eq.3) then
 +          y(1)=0.0D0
 +          y(2)=0.0D0
 +          else
 +        if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
  #ifdef OSF
            phii=phi(i)
-           icrc=0
-           call proc_proc(phii,icrc)
+ c          icrc=0
+ c          call proc_proc(phii,icrc)
            if (icrc.eq.1) phii=150.0
  #else
            phii=phi(i)
            y(1)=0.0D0
            y(2)=0.0D0
          endif
 -        if (i.lt.nres .and. itype(i).ne.ntyp1) then
 +        endif
 +        if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
  #ifdef OSF
            phii1=phi(i+1)
-           icrc=0
-           call proc_proc(phii1,icrc)
+ c          icrc=0
+ c          call proc_proc(phii1,icrc)
            if (icrc.eq.1) phii1=150.0
            phii1=pinorm(phii1)
            z(1)=cos(phii1)
        etheta=0.0D0
  c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
        do i=ithet_start,ithet_end
 +        if (i.le.2) cycle
 +        if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
 +     &  .or.itype(i).eq.ntyp1) cycle
+ c        if (itype(i-1).eq.ntyp1) cycle
 -        if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
 -     &(itype(i).eq.ntyp1)) cycle
          if (iabs(itype(i+1)).eq.20) iblock=2
          if (iabs(itype(i+1)).ne.20) iblock=1
          dethetai=0.0d0
Simple merge
Simple merge
Simple merge
@@@ -8,7 -8,7 +8,8 @@@
       & icheckgrad,minim,i2ndstr,refstr,pdbref,outpdb,outmol2,iprint,
       & overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb
       & ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file,
-      & constr_dist,gnorm_check,gradout,split_ene,symetr,AFMlog,
-      & selfguide
++     & selfguide,AFMlog,
+      & constr_dist,gnorm_check,gradout,split_ene,with_theta_constr,
+      & symetr
  C... minim = .true. means DO minimization.
  C... energy_dec = .true. means print energy decomposition matrix
Simple merge
@@@ -1,23 -1,23 +1,23 @@@
- FC = ifort
+ #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 -O0 -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
  
- CC = cc
- CFLAGS = -DLINUX -DPGI -c
- OPT =  -O3 -ip -w
- # -Mvect <---slows down
- #        -Minline=name:matmat2 <---false convergence
- LIBS = -Lxdrf -lxdrf
- #-DMOMENT
- #-DCO_BIAS
- #-DCRYST_TOR
- #-DDEBUG
+ #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
@@@ -42,43 -44,55 +42,56 @@@ object = unres.o arcos.o cartprint.o ch
          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 MREMD.o rattle.o gauss.o energy_split-sep.o \
 -      q_measure.o gnmr1.o ssMD.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 -DAMD64 -DUNRES -DISNAN \
-       -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC
- GAB: BIN = ../../../bin/unres/MD/unres-mult_ifort_single_GAB.exe
- GAB: ${object} xdrf/libxdrf.a
+ 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}
  
 -clean:
 -      /bin/rm *.o
 -
 -newconf.o: newconf.f
 -      ${FC} ${FFLAGS} ${CPPFLAGS} newconf.f
 -
 -bank.o: bank.F
 -      ${FC} ${FFLAGS} ${CPPFLAGS} bank.F
 -
 -diff12.o: diff12.f
 -      ${FC} ${FFLAGS} ${CPPFLAGS} diff12.f
 -
 -csa.o: csa.f
 -      ${FC} ${FFLAGS} ${CPPFLAGS} csa.f
 -
 -shift.o: shift.F
 -      ${FC} ${FFLAGS} ${CPPFLAGS} shift.F
 +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}
  
 -ran.o: ran.f
 -      ${FC} ${FFLAGS} ${CPPFLAGS} ran.f
 +xdrf/libxdrf.a:
 +      cd xdrf && make
  
 -together.o: together.F
 -      ${FC} ${FFLAGS} ${CPPFLAGS} together.F
 +clean:
 +      /bin/rm -f *.o && /bin/rm -f compinfo && cd xdrf && make clean
  
  test.o: test.F
        ${FC} ${FFLAGS} ${CPPFLAGS} test.F
Simple merge
        energia(17)=estr
        energia(20)=Uconst+Uconst_back
        energia(21)=esccor
 -C      energia(22)=eliptrans (the energy for lipid transfere implemented in lipid branch)
 -C      energia(23)= ... (energy for AFM, steered molecular dynamics)
 +      energia(22)=eliptran
 +      energia(23)=Eafmforce
+       energia(24)=ethetacnstr
  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"
@@@ -419,23 -396,22 +421,26 @@@ cMS$ATTRIBUTES C ::  proc_pro
        estr=energia(17)
        Uconst=energia(20)
        esccor=energia(21)
 +      eliptran=energia(22)
 +      Eafmforce=energia(23)
+       ethetacnstr=energia(24)
 -
  #ifdef SPLITELE
        etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
       & +wang*ebe+wtor*etors+wscloc*escloc
       & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
       & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
       & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
 -     & +wbond*estr+Uconst+wsccor*esccor+ethetacnstr
 +     & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce
++     & +ethetacnstr
  #else
        etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
       & +wang*ebe+wtor*etors+wscloc*escloc
       & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
       & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
       & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
 -     & +wbond*estr+Uconst+wsccor*esccor+ethetacnstr
 +     & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran
 +     & +Eafmforce
++     & +ethetacnstr
  #endif
        energia(0)=etot
  c detecting NaNQ
@@@ -1013,17 -976,17 +1018,18 @@@ C--------------------------------------
        estr=energia(17)
        Uconst=energia(20)
        esccor=energia(21)
 +      eliptran=energia(22)
 +      Eafmforce=energia(23) 
+       ethetacnstr=energia(24)
  #ifdef SPLITELE
        write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
       &  estr,wbond,ebe,wang,
       &  escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain,
       &  ecorr,wcorr,
       &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
--     &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
-      &  edihcnstr,ebr*nss,
-      &  Uconst,eliptran,wliptran,Eafmforce,etot
 -     &  edihcnstr,
 -     &  ethetacnstr,ebr*nss,
 -     &  Uconst,etot
++     &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
++     &  ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc,
++     &  etot
     10 format (/'Virtual-chain energies:'//
       & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
       & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
       & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
       & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
       & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
+      & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/
       & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
       & 'UCONST= ',1pE16.6,' (Constraint energy)'/ 
 +     & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
 +     & 'EAFM=  ',1pE16.6,' (atomic-force microscopy)'/
       & 'ETOT=  ',1pE16.6,' (total)')
 +
  #else
        write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,
       &  estr,wbond,ebe,wang,
       &  ecorr,wcorr,
       &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
       &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
-      &  ebr*nss,Uconst,eliptran,wliptran,Eafmforc,etot
 -     &  ethetacnstr,
 -     &  ebr*nss,Uconst,etot
++     &  ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc,
++     &  etot
     10 format (/'Virtual-chain energies:'//
       & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
       & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
       & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
       & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
       & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
+      & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/
       & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
       & 'UCONST=',1pE16.6,' (Constraint energy)'/ 
 +     & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
 +     & 'EAFM=  ',1pE16.6,' (atomic-force microscopy)'/
       & 'ETOT=  ',1pE16.6,' (total)')
  #endif
        return
        include 'COMMON.IOUNITS'
        include 'COMMON.CALC'
        include 'COMMON.CONTROL'
 +      include 'COMMON.SPLITELE'
        include 'COMMON.SBRIDGE'
        logical lprn
 -
 -c      write(iout,*) "Jestem w egb(evdw)"
 +      integer xshift,yshift,zshift
        evdw=0.0D0
  ccccc      energy_dec=.false.
 -c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 +C      print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=0.0D0
        lprn=.false.
  c     if (icall.eq.0) lprn=.false.
@@@ -5501,12 -4532,40 +5600,40 @@@ C Derivatives of the "mean" values in g
       &        E_theta,E_tc)
          endif
          etheta=etheta+ethetai
 -        if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
 -     &      'ebend',i,ethetai
 +        if (energy_dec) write (iout,'(a6,i5,0pf7.3,f7.3,i5)')
 +     &      'ebend',i,ethetai,theta(i),itype(i)
          if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
          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)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg)
        enddo
+       ethetacnstr=0.0d0
+ C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
+       do i=ithetaconstr_start,ithetaconstr_end
+         itheta=itheta_constr(i)
+         thetiii=theta(itheta)
+         difi=pinorm(thetiii-theta_constr0(i))
+         if (difi.gt.theta_drange(i)) then
+           difi=difi-theta_drange(i)
+           ethetacnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4
+           gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+      &    +for_thet_constr(i)*difi**3
+         else if (difi.lt.-drange(i)) then
+           difi=difi+drange(i)
+           ethetacnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4
+           gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+      &    +for_thet_constr(i)*difi**3
+         else
+           difi=0.0
+         endif
+        if (energy_dec) then
+         write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",
+      &    i,itheta,rad2deg*thetiii,
+      &    rad2deg*theta_constr0(i),  rad2deg*theta_drange(i),
+      &    rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
+      &    gloc(itheta+nphi-2,icg)
+         endif
+       enddo
  C Ufff.... We've done all this!!! 
        return
        end
        logical lprn /.false./, lprn1 /.false./
        etheta=0.0D0
        do i=ithet_start,ithet_end
 -        if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
 -     &(itype(i).eq.ntyp1)) cycle
 +c        print *,i,itype(i-1),itype(i),itype(i-2)
 +        if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
 +     &  .or.itype(i).eq.ntyp1) cycle
- C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
+ C        print *,i,theta(i)
          if (iabs(itype(i+1)).eq.20) iblock=2
          if (iabs(itype(i+1)).ne.20) iblock=1
          dethetai=0.0d0
            coskt(k)=dcos(k*theti2)
            sinkt(k)=dsin(k*theti2)
          enddo
+ C        print *,ethetai
 -
          if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
  #ifdef OSF
            phii=phi(i)
Simple merge
Simple merge
Simple merge
@@@ -71,8 -71,9 +71,9 @@@
          endif
        enddo
  #else
 -      read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
 +      read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok
        do i=1,ntyp
+       print *,i
          read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
       &   j=1,nbondterm(i)),msc(i),isc(i),restok(i)
          dsc(i) = vbldsc0(1,i)
          write (iout,*) "Coefficients of the cumulants"
        endif
        read (ifourier,*) nloctyp
 +
        do i=0,nloctyp-1
          read (ifourier,*,end=115,err=115)
 -        read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
 +        read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
 +#ifdef NEWCORR
 +        read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
 +        read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
 +        read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
 +        read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
 +        read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
 +#endif 
          if (lprint) then
          write (iout,*) 'Type',i
 -        write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 +        write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
          endif
 -        B1(1,i)  = b(3)
 -        B1(2,i)  = b(5)
 -        B1(1,-i) = b(3)
 -        B1(2,-i) = -b(5)
 +c        B1(1,i)  = b(3)
 +c        B1(2,i)  = b(5)
 +c        B1(1,-i) = b(3)
 +c        B1(2,-i) = -b(5)
  c        b1(1,i)=0.0d0
  c        b1(2,i)=0.0d0
 +c        B1tilde(1,i) = b(3)
 +c        B1tilde(2,i) =-b(5)
 +c        B1tilde(1,-i) =-b(3)
 +c        B1tilde(2,-i) =b(5)
 +c        b1tilde(1,i)=0.0d0
 +c        b1tilde(2,i)=0.0d0
 +c        B2(1,i)  = b(2)
 +c        B2(2,i)  = b(4)
 +c        B2(1,-i)  =b(2)
 +c        B2(2,-i)  =-b(4)
+         B1tilde(1,i) = b(3)
+         B1tilde(2,i) =-b(5)
+         B1tilde(1,-i) =-b(3)
+         B1tilde(2,-i) =b(5)
+         b1tilde(1,i)=0.0d0
+         b1tilde(2,i)=0.0d0
+         B2(1,i)  = b(2)
+         B2(2,i)  = b(4)
+         B2(1,-i)  =b(2)
+         B2(2,-i)  =-b(4)
  
  c        b2(1,i)=0.0d0
  c        b2(2,i)=0.0d0
@@@ -1108,20 -1092,12 +1119,20 @@@ C----------------------- LJK potential 
        goto 50
  C---------------------- GB or BP potential -----------------------------
     30 do i=1,ntyp
-        read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
+        read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
        enddo
 -      read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
 -      read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
 -      read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
 -      read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
 +      read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
 +      read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
 +      read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
 +      read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
 +C now we start reading lipid
 +      do i=1,ntyp
 +       read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
 +       print *,"WARNING!!"
 +       do j=1,ntyp
 +       epslip(i,j)=epslip(i,j)+0.05d0
 +       enddo
 +      enddo
  C For the GB potential convert sigma'**2 into chi'
        if (ipot.eq.4) then
        do i=1,ntyp
Simple merge
Simple merge
@@@ -1,30 -1,25 +1,31 @@@
        subroutine refsys(i2,i3,i4,e1,e2,e3,fail)
 +c This subroutine calculates unit vectors of a local reference system
 +c defined by atoms (i2), (i3), and (i4). The x axis is the axis from
+       implicit real*8 (a-h,o-z)
+       include 'DIMENSIONS'
+ c this subroutine calculates unity vectors of a local reference system
+ c defined by atoms (i2), (i3), and (i4). the x axis is the axis from
  c atom (i3) to atom (i2), and the xy plane is the plane defined by atoms
  c (i2), (i3), and (i4). z axis is directed according to the sign of the
- c vector product (i3)-(i2) and (i3)-(i4). Sets fail to .true. if atoms
+ c vector product (i3)-(i2) and (i3)-(i4). sets fail to .true. if atoms
  c (i2) and (i3) or (i3) and (i4) coincide or atoms (i2), (i3), and (i4)
- c form a linear fragment. Returns vectors e1, e2, and e3.
-       implicit real*8 (a-h,o-z)
-       include 'DIMENSIONS'
+ c form a linear fragment. returns vectors e1, e2, and e3.
        logical fail
        double precision e1(3),e2(3),e3(3)
        double precision u(3),z(3)
        include 'COMMON.IOUNITS'
-       include 'COMMON.CHAIN'
 -      include "COMMON.CHAIN"
 -      data coinc /1.0d-13/,align /1.0d-13/
 +      double precision coinc/1.0D-13/,align /1.0D-13/
 +c      print *,'just initialize'
        fail=.false.
 -      s1=0.0d0
 -      s2=0.0d0
 +c      print *,fail
 +      s1=0.0
 +      s2=0.0
 +      print *,s1,s2
        do 1 i=1,3
 +      print *, i2,i3,i4
        zi=c(i,i2)-c(i,i3)
        ui=c(i,i4)-c(i,i3)
 +      print *,zi,ui
        s1=s1+zi*zi
        s2=s2+ui*ui
        z(i)=zi
      2 if (s2.gt.coinc) goto 4
        write(iout,1000) i3,i4,i1
        fail=.true.
-       do 5 i=1,3
-     5 c(i,i1)=0.0D0
        return
 +      print *,'two if pass'
      4 s1=1.0/s1
        s2=1.0/s2
        v1=z(2)*u(3)-z(3)*u(2)
        e2(1)=e1(3)*e3(2)-e1(2)*e3(3)
        e2(2)=e1(1)*e3(3)-e1(3)*e3(1)
        e2(3)=e1(2)*e3(1)-e1(1)*e3(2)
-       print *,'just before leave'
 - 1000 format (/1x,' * * * error - atoms',i4,' and',i4,' coincide.')
 - 1010 format (/1x,' * * * error - atoms',2(i4,2h, ),i4,' form a linear')
 + 1000 format (/1x,' * * * Error - atoms',i4,' and',i4,' coincide.',
 +     1 'coordinates of atom',i4,' are set to zero.')
 + 1010 format (/1x,' * * * Error - atoms',2(i4,2h, ),i4,' form a linear',
 +     1 ' fragment. coordinates of atom',i4,' are set to zero.')
        return
        end
Simple merge
Simple merge
Simple merge
@@@ -210,10 -210,9 +210,11 @@@ C Kozlowska et al., J. Phys.: Condens. 
  C
        read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
       &  ntheterm3,nsingle,ndouble
+ C      print *, "tu"
        nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 +      write(iout,*) "I am here",ntyp1
        read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
 +C      write(iout,*) "I am herew"
        do i=-ntyp1,-1
          ithetyp(i)=-ithetyp(-i)
        enddo
@@@ -660,19 -644,15 +662,19 @@@ c      &v2(k,-i,-j,iblock),v2(k,i,j,ibl
  C
  C 6/23/01 Read parameters for double torsionals
  C
 -      do i=1,ntortyp
 -        do j=1,ntortyp
 -          do k=1,ntortyp
 +C      do i=1,ntortyp
 +C        do j=1,ntortyp
 +C          do k=1,ntortyp
 +      do iblock=1,2
 +      do i=0,ntortyp-1
 +        do j=-ntortyp+1,ntortyp-1
 +          do k=-ntortyp+1,ntortyp-1
              read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
- c              write (iout,*) "OK onelett",
- c     &         i,j,k,t1,t2,t3
+               write (iout,*) "OK onelett",
+      &         i,j,k,t1,t2,t3
  
-             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) 
-      &        .or. t3.ne.toronelet(k)) then
+             if (t1.ne.onelett(i) .or. t2.ne.onelett(j) 
+      &        .or. t3.ne.onelett(k)) then
                write (iout,*) "Error in double torsional parameter file",
       &         i,j,k,t1,t2,t3
  #ifdef MPI
Simple merge
Simple merge
@@@ -70,8 -67,8 +70,9 @@@ cd    print *,'EHPB exitted succesfully
  C
  C Calculate the virtual-bond-angle energy.
  C
-       call ebend(ebe)
 +C      print *,'Bend energy finished.'
+       call ebend(ebe,ethetacnstr)
+ cd    print *,'Bend energy finished.'
  C
  C Calculate the SC local energy.
  C
@@@ -3390,16 -3315,12 +3504,16 @@@ C Zero the energy function and its deri
            ichir21=isign(1,itype(i))
            ichir22=isign(1,itype(i))
           endif
 +         if (i.eq.3) then
 +          y(1)=0.0D0
 +          y(2)=0.0D0
 +          else
  
 -        if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
 +        if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
  #ifdef OSF
            phii=phi(i)
-           icrc=0
-           call proc_proc(phii,icrc)
+ c          icrc=0
+ c          call proc_proc(phii,icrc)
            if (icrc.eq.1) phii=150.0
  #else
            phii=phi(i)
            y(1)=0.0D0
            y(2)=0.0D0
          endif
 -        if (i.lt.nres .and. itype(i).ne.ntyp1) then
 +        endif
 +        if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
  #ifdef OSF
            phii1=phi(i+1)
-           icrc=0
-           call proc_proc(phii1,icrc)
+ c          icrc=0
+ c          call proc_proc(phii1,icrc)
            if (icrc.eq.1) phii1=150.0
            phii1=pinorm(phii1)
            z(1)=cos(phii1)
Simple merge
Simple merge
@@@ -17,7 -17,7 +17,8 @@@
        include "COMMON.FREE"
        include "COMMON.CONTROL"
        include "COMMON.ENERGIES"
 +      include "COMMON.SPLITELE"
+       include "COMMON.SBRIDGE"
        character*800 controlcard
        integer i,j,k,ii,n_ene_found
        integer ind,itype1,itype2,itypf,itypsc,itypp