Fixed the addition of dummy chain ends to structures read from UNRES PDB files
authorAdam Liwo <jal47@matrix.chem.cornell.edu>
Sun, 17 Feb 2013 03:34:55 +0000 (22:34 -0500)
committerAdam Liwo <jal47@matrix.chem.cornell.edu>
Sun, 17 Feb 2013 04:04:41 +0000 (23:04 -0500)
Fixed the bug in arcos.f (returned alpha=0 for cos(alpha)<=-1)

Conflicts:

bin/unres/MD/unres_ifort_MPICH_E0LL2Y.exe
source/unres/src_MD-M/Makefile

30 files changed:
bin/unres/CSA/unres_csa_ifort_mpich-1.2.7p1.exe
bin/unres/CSA/unres_dfa_csa-Yi.exe
bin/unres/MD/unres_Tc_procor_oldparm_em64-D-symetr.exe
bin/unres/MINIM/unres_ifort_MIN_single_E0LL2Y.exe
bin/wham/wham_ifort_MPICH_E0LL2Y-DEBUG.exe
source/unres/src_CSA/Makefile_GAB [new file with mode: 0644]
source/unres/src_CSA/arcos.f
source/unres/src_CSA/readpdb.F
source/unres/src_CSA/refsys.f [new file with mode: 0644]
source/unres/src_CSA_DiL/arcos.f
source/unres/src_CSA_DiL/readpdb.F
source/unres/src_MD-M/COMMON.REFSYS [deleted file]
source/unres/src_MD-M/Makefile
source/unres/src_MD-M/arcos.f
source/unres/src_MD-M/bond_move.f
source/unres/src_MD-M/readpdb.F
source/unres/src_MD-M/refsys.f
source/unres/src_MD/COMMON.REFSYS [deleted file]
source/unres/src_MD/arcos.f
source/unres/src_MD/bond_move.f
source/unres/src_MD/readpdb.F
source/unres/src_MD/refsys.f
source/unres/src_MIN/Makefile [new symlink]
source/unres/src_MIN/Makefile_ifort_single
source/unres/src_MIN/arcos.f
source/unres/src_MIN/readpdb.F
source/unres/src_MIN/refsys.f [new file with mode: 0644]
source/wham/src/enecalc1.F
source/wham/src/molread_zs.F
source/wham/src/readrtns.F

index 6710349..fca8c0a 100755 (executable)
Binary files a/bin/unres/CSA/unres_csa_ifort_mpich-1.2.7p1.exe and b/bin/unres/CSA/unres_csa_ifort_mpich-1.2.7p1.exe differ
index 528b3aa..2efc549 100755 (executable)
Binary files a/bin/unres/CSA/unres_dfa_csa-Yi.exe and b/bin/unres/CSA/unres_dfa_csa-Yi.exe differ
index 5d0105f..932579c 100755 (executable)
Binary files a/bin/unres/MD/unres_Tc_procor_oldparm_em64-D-symetr.exe and b/bin/unres/MD/unres_Tc_procor_oldparm_em64-D-symetr.exe differ
index cb1a624..78e6748 100755 (executable)
Binary files a/bin/unres/MINIM/unres_ifort_MIN_single_E0LL2Y.exe and b/bin/unres/MINIM/unres_ifort_MIN_single_E0LL2Y.exe differ
index a9af08d..b58e33d 100755 (executable)
Binary files a/bin/wham/wham_ifort_MPICH_E0LL2Y-DEBUG.exe and b/bin/wham/wham_ifort_MPICH_E0LL2Y-DEBUG.exe differ
diff --git a/source/unres/src_CSA/Makefile_GAB b/source/unres/src_CSA/Makefile_GAB
new file mode 100644 (file)
index 0000000..44d2419
--- /dev/null
@@ -0,0 +1,99 @@
+CPPFLAGS = -DPROCOR -DLINUX -DPGI -DISNAN -DMP -DMPI -DUNRES \
+           -DSPLITELE -DAMD64 -DLANG0 -DPROCOR \
+           -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DPROCOR
+#           -DTSCSC
+#-DTIMING \
+# -DCRYST_BOND -DCRYST_THETA -DCRYST_SC 
+# -DMOMENT
+#-DPARVEC 
+#-DPARINT -DPARINTDER  
+
+INSTALL_DIR = /users/software/mpich-1.2.7p1_intel-10.1_em64_ssh
+
+FC= ifort 
+
+OPT =  -O3 -ip -w 
+
+FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include
+FFLAGS1 = -c -w -g -O0 -d2 -CA -CB -I$(INSTALL_DIR)/include
+FFLAGS2 = -c -w -g -O0 -I$(INSTALL_DIR)/include
+FFLAGS3 = -c -w -O3 -mp
+FFLAGSE = -c -w -O3 -ipo -ipo_obj  -opt_report -I$(INSTALL_DIR)/include
+
+
+BIN = ../../../bin/unres/CSA/unres_csa_ifort_mpich-1.2.7p1.exe
+LIBS =  -lpthread -L$(INSTALL_DIR)/lib -lmpich
+
+ARCH = LINUX
+PP = /lib/cpp -P
+
+
+all: unres
+
+.SUFFIXES: .F
+.F.o:
+       ${FC} ${FFLAGS} ${CPPFLAGS} $*.F
+
+
+object = unres_csa.o arcos.o cartprint.o chainbuild.o initialize_p.o \
+        matmult.o readrtns_csa.o parmread.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 \
+       gradient_p.o minimize_p.o sumsld.o \
+        cored.o rmdd.o geomout_min.o readpdb.o \
+        intcartderiv.o \
+        MP.o printmat.o convert.o int_to_cart.o \
+       dfa.o \
+        together.o csa.o minim_jlee.o shift.o diff12.o bank.o newconf.o ran.o \
+        indexx.o prng_32.o contact.o gen_rand_conf.o \
+        sc_move.o test.o local_move.o rmsd.o fitsq.o elecont.o djacob.o \
+        distfit.o banach.o TMscore_subroutine.o minim_mult.o refsys.o
+
+unres: ${object} 
+       cc -o compinfo compinfo.c
+       ./compinfo | true
+       ${FC} ${FFLAGS} cinfo.f
+       ${FC} ${OPT} ${object} cinfo.o ${LIBS}  -o ${BIN}
+
+
+clean:
+       /bin/rm *.o *.il
+
+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} ${FFLAGS3} ${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
+
+dfa.o : dfa.F
+       ${FC} -mp  ${FFLAGS3} ${CPPFLAGS} dfa.F
+
+
+
index 69810ea..f054118 100644 (file)
@@ -2,7 +2,7 @@
       implicit real*8 (a-h,o-z)
       include 'COMMON.GEO'
       IF (DABS(X).LT.1.0D0) GOTO 1
-      ARCOS=0.5D0*(PI+DSIGN(1.0D0,X)*PI)
+      ARCOS=PIPOL*(1.0d0-DSIGN(1.0D0,X))
       RETURN
     1 ARCOS=DACOS(X)
       RETURN
index 69ec3e1..eb4ba3f 100644 (file)
@@ -16,6 +16,8 @@ C geometry.
       character*3 seq,atom,res
       character*80 card
       dimension sccor(3,20)
+      double precision e1(3),e2(3),e3(3)
+      logical fail
       integer rescode
       ibeg=1
       lsecondary=.false.
@@ -104,9 +106,16 @@ C Calculate the CM of the last side chain.
         nres=nres+1
         itype(nres)=21
         if (unres_pdb) then
-          c(1,nres)=c(1,nres-1)+3.8d0
-          c(2,nres)=c(2,nres-1)
-          c(3,nres)=c(3,nres-1)
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+          call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+          if (fail) then
+            e2(1)=0.0d0
+            e2(2)=1.0d0
+            e2(3)=0.0d0
+          endif
+          do j=1,3
+            c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
+          enddo
         else
         do j=1,3
           dcj=c(j,nres-2)-c(j,nres-3)
@@ -128,9 +137,16 @@ C Calculate the CM of the last side chain.
         nsup=nsup-1
         nstart_sup=2
         if (unres_pdb) then
-          c(1,1)=c(1,2)-3.8d0
-          c(2,1)=c(2,2)
-          c(3,1)=c(3,2)
+C 2/15/2013 by Adam: corrected insertion of the first dummy residue
+          call refsys(2,3,4,e1,e2,e3,fail)
+          if (fail) then
+            e2(1)=0.0d0
+            e2(2)=1.0d0
+            e2(3)=0.0d0
+          endif
+          do j=1,3
+            c(j,1)=c(j,2)-3.8d0*e2(j)
+          enddo
         else
         do j=1,3
           dcj=c(j,4)-c(j,3)
@@ -258,8 +274,8 @@ c      endif
           enddo
           iti=itype(i)
           di=dist(i,nres+i)
-C 10/03/12 Adam: Correction for zero SC-SC bond length
-          if (itype(i).ne.10 .and. itype(i).ne.21. and. di.eq.0.0d0)
+C 9/29/12 Adam: Correction for zero SC-SC bond length 
+          if (itype(i).ne.10 .and. itype(i).ne.21. and. di.eq.0.0d0) 
      &     di=dsc(itype(i))
           vbld(i+nres)=di
           if (itype(i).ne.10) then
@@ -410,4 +426,3 @@ c       print *,vbld(i+1),vbld(i+1+nres)
       enddo
       return
       end
-      
diff --git a/source/unres/src_CSA/refsys.f b/source/unres/src_CSA/refsys.f
new file mode 100644 (file)
index 0000000..b57c201
--- /dev/null
@@ -0,0 +1,60 @@
+      subroutine refsys(i2,i3,i4,e1,e2,e3,fail)
+      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 (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.
+      logical fail
+      double precision e1(3),e2(3),e3(3)
+      double precision u(3),z(3)
+      include 'COMMON.IOUNITS'
+      include "COMMON.CHAIN"
+      data coinc /1.0d-13/,align /1.0d-13/
+      fail=.false.
+      s1=0.0d0
+      s2=0.0d0
+      do 1 i=1,3
+      zi=c(i,i2)-c(i,i3)
+      ui=c(i,i4)-c(i,i3)
+      s1=s1+zi*zi
+      s2=s2+ui*ui
+      z(i)=zi
+    1 u(i)=ui
+      s1=sqrt(s1)
+      s2=sqrt(s2)
+      if (s1.gt.coinc) goto 2
+      write (iout,1000) i2,i3,i1
+      fail=.true.
+      return
+    2 if (s2.gt.coinc) goto 4
+      write(iout,1000) i3,i4,i1
+      fail=.true.
+      return
+    4 s1=1.0/s1
+      s2=1.0/s2
+      v1=z(2)*u(3)-z(3)*u(2)
+      v2=z(3)*u(1)-z(1)*u(3)
+      v3=z(1)*u(2)-z(2)*u(1)
+      anorm=sqrt(v1*v1+v2*v2+v3*v3)
+      if (anorm.gt.align) goto 6
+      write (iout,1010) i2,i3,i4,i1
+      fail=.true.
+      return
+    6 anorm=1.0/anorm
+      e3(1)=v1*anorm
+      e3(2)=v2*anorm
+      e3(3)=v3*anorm
+      e1(1)=z(1)*s1
+      e1(2)=z(2)*s1
+      e1(3)=z(3)*s1
+      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)
+ 1000 format (/1x,' * * * error - atoms',i4,' and',i4,' coincide.')
+ 1010 format (/1x,' * * * error - atoms',2(i4,2h, ),i4,' form a linear')
+      return
+      end
index 69810ea..f054118 100644 (file)
@@ -2,7 +2,7 @@
       implicit real*8 (a-h,o-z)
       include 'COMMON.GEO'
       IF (DABS(X).LT.1.0D0) GOTO 1
-      ARCOS=0.5D0*(PI+DSIGN(1.0D0,X)*PI)
+      ARCOS=PIPOL*(1.0d0-DSIGN(1.0D0,X))
       RETURN
     1 ARCOS=DACOS(X)
       RETURN
index 69ec3e1..7810b24 100644 (file)
@@ -16,6 +16,8 @@ C geometry.
       character*3 seq,atom,res
       character*80 card
       dimension sccor(3,20)
+      double precision e1(3),e2(3),e3(3)
+      logical fail
       integer rescode
       ibeg=1
       lsecondary=.false.
@@ -104,9 +106,16 @@ C Calculate the CM of the last side chain.
         nres=nres+1
         itype(nres)=21
         if (unres_pdb) then
-          c(1,nres)=c(1,nres-1)+3.8d0
-          c(2,nres)=c(2,nres-1)
-          c(3,nres)=c(3,nres-1)
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+          call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+          if (fail) then
+            e2(1)=0.0d0
+            e2(2)=1.0d0
+            e2(3)=0.0d0
+          endif
+          do j=1,3
+            c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
+          enddo
         else
         do j=1,3
           dcj=c(j,nres-2)-c(j,nres-3)
@@ -128,9 +137,16 @@ C Calculate the CM of the last side chain.
         nsup=nsup-1
         nstart_sup=2
         if (unres_pdb) then
-          c(1,1)=c(1,2)-3.8d0
-          c(2,1)=c(2,2)
-          c(3,1)=c(3,2)
+C 2/15/2013 by Adam: corrected insertion of the first dummy residue
+          call refsys(2,3,4,e1,e2,e3,fail)
+          if (fail) then
+            e2(1)=0.0d0
+            e2(2)=1.0d0
+            e2(3)=0.0d0
+          endif
+          do j=1,3
+            c(j,1)=c(j,2)-3.8d0*e2(j)
+          enddo
         else
         do j=1,3
           dcj=c(j,4)-c(j,3)
@@ -258,8 +274,8 @@ c      endif
           enddo
           iti=itype(i)
           di=dist(i,nres+i)
-C 10/03/12 Adam: Correction for zero SC-SC bond length
-          if (itype(i).ne.10 .and. itype(i).ne.21. and. di.eq.0.0d0)
+C 9/29/12 Adam: Correction for zero SC-SC bond length 
+          if (itype(i).ne.10 .and. itype(i).ne.21. and. di.eq.0.0d0) 
      &     di=dsc(itype(i))
           vbld(i+nres)=di
           if (itype(i).ne.10) then
@@ -410,4 +426,64 @@ c       print *,vbld(i+1),vbld(i+1+nres)
       enddo
       return
       end
+C-----------------------------------------------------------------------------
+      subroutine refsys(i2,i3,i4,e1,e2,e3,fail)
+      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 (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.
+      logical fail
+      double precision e1(3),e2(3),e3(3)
+      double precision u(3),z(3)
+      include "COMMON.CHAIN"
+      data coinc /1.0d-13/,align /1.0d-13/
+      fail=.false.
+      s1=0.0d0
+      s2=0.0d0
+      do 1 i=1,3
+      zi=c(i,i2)-c(i,i3)
+      ui=c(i,i4)-c(i,i3)
+      s1=s1+zi*zi
+      s2=s2+ui*ui
+      z(i)=zi
+    1 u(i)=ui
+      s1=sqrt(s1)
+      s2=sqrt(s2)
+      if (s1.gt.coinc) goto 2
+      write (iout,1000) i2,i3,i1
+      fail=.true.
+      return
+    2 if (s2.gt.coinc) goto 4
+      write(iout,1000) i3,i4,i1
+      fail=.true.
+      return
+    4 s1=1.0/s1
+      s2=1.0/s2
+      v1=z(2)*u(3)-z(3)*u(2)
+      v2=z(3)*u(1)-z(1)*u(3)
+      v3=z(1)*u(2)-z(2)*u(1)
+      anorm=sqrt(v1*v1+v2*v2+v3*v3)
+      if (anorm.gt.align) goto 6
+      write (iout,1010) i2,i3,i4,i1
+      fail=.true.
+      return
+    6 anorm=1.0/anorm
+      e3(1)=v1*anorm
+      e3(2)=v2*anorm
+      e3(3)=v3*anorm
+      e1(1)=z(1)*s1
+      e1(2)=z(2)*s1
+      e1(3)=z(3)*s1
+      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)
+ 1000 format (/1x,' * * * error - atoms',i4,' and',i4,' coincide.')
+ 1010 format (/1x,' * * * error - atoms',2(i4,2h, ),i4,' form a linear')
+      return
+      end
       
diff --git a/source/unres/src_MD-M/COMMON.REFSYS b/source/unres/src_MD-M/COMMON.REFSYS
deleted file mode 100644 (file)
index 9eaa3c3..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-      double precision e1,e2,e3,u,z,s1,s2
-      integer i1,i2,i3,i4
-      common /refer/ e1(3),e2(3),e3(3),u(3),z(3),s1,s2,i1,i2,i3,i4
index c850dec..8db4abe 100644 (file)
@@ -23,8 +23,8 @@ FFLAGS1 = -c -w -g -d2 -CA -CB -I$(INSTALL_DIR)/include
 FFLAGS2 = -c -w -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_oldparm_em64-D-symetr.exe
-#LIBS = -L$(INSTALL_DIR)/lib_pgi -lmpich xdrf/libxdrf.a
+BIN = ../../../bin/unres/MD/unres_Tc_procor_oldparm_em64-D-symetr.exe
+#LIBS = -L$(INSTALL_DIR)/lib_pgi -lmpich ../../lib/xdrf/libxdrf.a
 #LIBS = -L$(INSTALL_DIR)/lib_ifort -lmpich xdrf/libxdrf.a
 LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a -g -d2 -CA -CB
 
@@ -60,7 +60,8 @@ unres: ${object} proc_proc.o
        cc -o compinfo compinfo.c 
        ./compinfo | true
        ${FC} ${FFLAGS} cinfo.f
-       ${FC} ${OPT} -Wl,-M ${object} proc_proc.o cinfo.o ${LIBS}  -o ${BIN}
+       ${FC} ${OPT} ${object} proc_proc.o cinfo.o ${LIBS}  -o ${BIN}
+#${FC} ${OPT} -Wl,-M ${object} proc_proc.o cinfo.o ${LIBS}  -o ${BIN}
 
 
 clean:
@@ -128,3 +129,12 @@ lagrangian_lesyng.o : lagrangian_lesyng.F
 
 proc_proc.o: proc_proc.c
        ${CC} ${CFLAGS} proc_proc.c
+
+add.o: add.f
+       ${FC} ${FFLAGS2} add.f
+
+blas.o: blas.f
+       ${FC} ${FFLAGS2} blas.f
+
+eigen.o: eigen.f
+       ${FC} ${FFLAGS2} eigen.f
index 69810ea..f054118 100644 (file)
@@ -2,7 +2,7 @@
       implicit real*8 (a-h,o-z)
       include 'COMMON.GEO'
       IF (DABS(X).LT.1.0D0) GOTO 1
-      ARCOS=0.5D0*(PI+DSIGN(1.0D0,X)*PI)
+      ARCOS=PIPOL*(1.0d0-DSIGN(1.0D0,X))
       RETURN
     1 ARCOS=DACOS(X)
       RETURN
index 4c0761a..4843f60 100644 (file)
@@ -8,10 +8,9 @@ C Move NBOND fragment starting from the CA(nstart) by angle PSI.
       include 'COMMON.GEO'
       include 'COMMON.CHAIN'
       include 'COMMON.VAR'
-      include 'COMMON.REFSYS'
       include 'COMMON.IOUNITS'
       include 'COMMON.MCM'
-      dimension x(3),e(3,3),rot(3,3),trans(3,3)
+      dimension x(3),e(3,3),e1(3),e2(3),e3(3),rot(3,3),trans(3,3)
       error=.false.
       nend=nstart+nbond
       if (print_mc.gt.2) then
@@ -32,7 +31,7 @@ C Generate the reference system.
       i2=nend
       i3=nstart
       i4=nstart+1
-      call refsys(error) 
+      call refsys(i2,i3,i4,e1,e2,e3,error) 
 C Return, if couldn't define the reference system.
       if (error) return
 C Compute the transformation matrix.
index edd478e..8b6f331 100644 (file)
@@ -16,7 +16,9 @@ C geometry.
       character*3 seq,atom,res
       character*80 card
       dimension sccor(3,20)
+      double precision e1(3),e2(3),e3(3)
       integer rescode
+      logical fail
       ibeg=1
       lsecondary=.false.
       nhfrag=0
@@ -140,9 +142,16 @@ C Calculate the CM of the last side chain.
         nres=nres+1
         itype(nres)=21
         if (unres_pdb) then
-          c(1,nres)=c(1,nres-1)+3.8d0
-          c(2,nres)=c(2,nres-1)
-          c(3,nres)=c(3,nres-1)
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+          call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+          if (fail) then
+            e2(1)=0.0d0
+            e2(2)=1.0d0
+            e2(3)=0.0d0
+          endif
+          do j=1,3
+            c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
+          enddo
         else
         do j=1,3
           dcj=c(j,nres-2)-c(j,nres-3)
@@ -164,9 +173,16 @@ C Calculate the CM of the last side chain.
         nsup=nsup-1
         nstart_sup=2
         if (unres_pdb) then
-          c(1,1)=c(1,2)-3.8d0
-          c(2,1)=c(2,2)
-          c(3,1)=c(3,2)
+C 2/15/2013 by Adam: corrected insertion of the first dummy residue
+          call refsys(2,3,4,e1,e2,e3,fail)
+          if (fail) then
+            e2(1)=0.0d0
+            e2(2)=1.0d0
+            e2(3)=0.0d0
+          endif
+          do j=1,3
+            c(j,1)=c(j,2)-3.8d0*e2(j)
+          enddo
         else
         do j=1,3
           dcj=c(j,4)-c(j,3)
index ec620df..86b0524 100644 (file)
@@ -9,10 +9,11 @@ c form a linear fragment. Returns vectors e1, e2, and e3.
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       logical fail
+      double precision e1(3),e2(3),e3(3)
+      double precision u(3),z(3)
       include 'COMMON.IOUNITS'
       include 'COMMON.CHAIN'
-      include 'COMMON.REFSYS'
-      double precision coinc/1.0D-4/,align /1.0D-7/
+      double precision coinc/1.0D-13/,align /1.0D-13/
       fail=.false.
       s1=0.0
       s2=0.0
diff --git a/source/unres/src_MD/COMMON.REFSYS b/source/unres/src_MD/COMMON.REFSYS
deleted file mode 100644 (file)
index 9eaa3c3..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-      double precision e1,e2,e3,u,z,s1,s2
-      integer i1,i2,i3,i4
-      common /refer/ e1(3),e2(3),e3(3),u(3),z(3),s1,s2,i1,i2,i3,i4
index 69810ea..f054118 100644 (file)
@@ -2,7 +2,7 @@
       implicit real*8 (a-h,o-z)
       include 'COMMON.GEO'
       IF (DABS(X).LT.1.0D0) GOTO 1
-      ARCOS=0.5D0*(PI+DSIGN(1.0D0,X)*PI)
+      ARCOS=PIPOL*(1.0d0-DSIGN(1.0D0,X))
       RETURN
     1 ARCOS=DACOS(X)
       RETURN
index 4c0761a..4843f60 100644 (file)
@@ -8,10 +8,9 @@ C Move NBOND fragment starting from the CA(nstart) by angle PSI.
       include 'COMMON.GEO'
       include 'COMMON.CHAIN'
       include 'COMMON.VAR'
-      include 'COMMON.REFSYS'
       include 'COMMON.IOUNITS'
       include 'COMMON.MCM'
-      dimension x(3),e(3,3),rot(3,3),trans(3,3)
+      dimension x(3),e(3,3),e1(3),e2(3),e3(3),rot(3,3),trans(3,3)
       error=.false.
       nend=nstart+nbond
       if (print_mc.gt.2) then
@@ -32,7 +31,7 @@ C Generate the reference system.
       i2=nend
       i3=nstart
       i4=nstart+1
-      call refsys(error) 
+      call refsys(i2,i3,i4,e1,e2,e3,error) 
 C Return, if couldn't define the reference system.
       if (error) return
 C Compute the transformation matrix.
index 97e9aa8..48e0abd 100644 (file)
@@ -16,6 +16,8 @@ C geometry.
       character*3 seq,atom,res
       character*80 card
       dimension sccor(3,20)
+      double precision e1(3),e2(3),e3(3)
+      logical fail
       integer rescode
       ibeg=1
       lsecondary=.false.
@@ -104,9 +106,16 @@ C Calculate the CM of the last side chain.
         nres=nres+1
         itype(nres)=21
         if (unres_pdb) then
-          c(1,nres)=c(1,nres-1)+3.8d0
-          c(2,nres)=c(2,nres-1)
-          c(3,nres)=c(3,nres-1)
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+          call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+          if (fail) then
+            e2(1)=0.0d0
+            e2(2)=1.0d0
+            e2(3)=0.0d0
+          endif
+          do j=1,3
+            c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
+          enddo
         else
         do j=1,3
           dcj=c(j,nres-2)-c(j,nres-3)
@@ -128,9 +137,16 @@ C Calculate the CM of the last side chain.
         nsup=nsup-1
         nstart_sup=2
         if (unres_pdb) then
-          c(1,1)=c(1,2)-3.8d0
-          c(2,1)=c(2,2)
-          c(3,1)=c(3,2)
+C 2/15/2013 by Adam: corrected insertion of the first dummy residue
+          call refsys(2,3,4,e1,e2,e3,fail)
+          if (fail) then
+            e2(1)=0.0d0
+            e2(2)=1.0d0
+            e2(3)=0.0d0
+          endif
+          do j=1,3
+            c(j,1)=c(j,2)-3.8d0*e2(j)
+          enddo
         else
         do j=1,3
           dcj=c(j,4)-c(j,3)
@@ -414,4 +430,3 @@ c       print *,vbld(i+1),vbld(i+1+nres)
       enddo
       return
       end
-      
index ec620df..b57c201 100644 (file)
@@ -1,21 +1,22 @@
-      subroutine refsys(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
+      subroutine refsys(i2,i3,i4,e1,e2,e3,fail)
+      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.REFSYS'
-      double precision coinc/1.0D-4/,align /1.0D-7/
+      include "COMMON.CHAIN"
+      data coinc /1.0d-13/,align /1.0d-13/
       fail=.false.
-      s1=0.0
-      s2=0.0
+      s1=0.0d0
+      s2=0.0d0
       do 1 i=1,3
       zi=c(i,i2)-c(i,i3)
       ui=c(i,i4)-c(i,i3)
@@ -28,28 +29,22 @@ c form a linear fragment. Returns vectors e1, e2, and e3.
       if (s1.gt.coinc) goto 2
       write (iout,1000) i2,i3,i1
       fail=.true.
-c     do 3 i=1,3
-c   3 c(i,i1)=0.0D0
       return
     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
     4 s1=1.0/s1
       s2=1.0/s2
       v1=z(2)*u(3)-z(3)*u(2)
       v2=z(3)*u(1)-z(1)*u(3)
       v3=z(1)*u(2)-z(2)*u(1)
-      anorm=dsqrt(v1*v1+v2*v2+v3*v3)
+      anorm=sqrt(v1*v1+v2*v2+v3*v3)
       if (anorm.gt.align) goto 6
       write (iout,1010) i2,i3,i4,i1
       fail=.true.
-c     do 7 i=1,3
-c   7 c(i,i1)=0.0D0
       return
-    6 anorm=1.0D0/anorm
+    6 anorm=1.0/anorm
       e3(1)=v1*anorm
       e3(2)=v2*anorm
       e3(3)=v3*anorm
@@ -59,9 +54,7 @@ c   7 c(i,i1)=0.0D0
       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)
- 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.')
+ 1000 format (/1x,' * * * error - atoms',i4,' and',i4,' coincide.')
+ 1010 format (/1x,' * * * error - atoms',2(i4,2h, ),i4,' form a linear')
       return
       end
diff --git a/source/unres/src_MIN/Makefile b/source/unres/src_MIN/Makefile
new file mode 120000 (symlink)
index 0000000..cbb1cd3
--- /dev/null
@@ -0,0 +1 @@
+Makefile_ifort_single
\ No newline at end of file
index d21cc98..1e5d224 100644 (file)
@@ -30,7 +30,7 @@ object = unres_min.o arcos.o cartprint.o chainbuild.o initialize_p.o \
         cored.o rmdd.o geomout_min.o readpdb.o \
         intcartderiv.o \
         MP.o printmat.o convert.o int_to_cart.o \
-       djacob.o gen_rand_conf.o sc_move.o
+       djacob.o gen_rand_conf.o sc_move.o refsys.o
 
 GAB: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNRES -DISNAN \
        -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC
index 69810ea..f054118 100644 (file)
@@ -2,7 +2,7 @@
       implicit real*8 (a-h,o-z)
       include 'COMMON.GEO'
       IF (DABS(X).LT.1.0D0) GOTO 1
-      ARCOS=0.5D0*(PI+DSIGN(1.0D0,X)*PI)
+      ARCOS=PIPOL*(1.0d0-DSIGN(1.0D0,X))
       RETURN
     1 ARCOS=DACOS(X)
       RETURN
index 1d79870..eb4ba3f 100644 (file)
@@ -16,6 +16,8 @@ C geometry.
       character*3 seq,atom,res
       character*80 card
       dimension sccor(3,20)
+      double precision e1(3),e2(3),e3(3)
+      logical fail
       integer rescode
       ibeg=1
       lsecondary=.false.
@@ -104,9 +106,16 @@ C Calculate the CM of the last side chain.
         nres=nres+1
         itype(nres)=21
         if (unres_pdb) then
-          c(1,nres)=c(1,nres-1)+3.8d0
-          c(2,nres)=c(2,nres-1)
-          c(3,nres)=c(3,nres-1)
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+          call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+          if (fail) then
+            e2(1)=0.0d0
+            e2(2)=1.0d0
+            e2(3)=0.0d0
+          endif
+          do j=1,3
+            c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
+          enddo
         else
         do j=1,3
           dcj=c(j,nres-2)-c(j,nres-3)
@@ -128,9 +137,16 @@ C Calculate the CM of the last side chain.
         nsup=nsup-1
         nstart_sup=2
         if (unres_pdb) then
-          c(1,1)=c(1,2)-3.8d0
-          c(2,1)=c(2,2)
-          c(3,1)=c(3,2)
+C 2/15/2013 by Adam: corrected insertion of the first dummy residue
+          call refsys(2,3,4,e1,e2,e3,fail)
+          if (fail) then
+            e2(1)=0.0d0
+            e2(2)=1.0d0
+            e2(3)=0.0d0
+          endif
+          do j=1,3
+            c(j,1)=c(j,2)-3.8d0*e2(j)
+          enddo
         else
         do j=1,3
           dcj=c(j,4)-c(j,3)
@@ -410,4 +426,3 @@ c       print *,vbld(i+1),vbld(i+1+nres)
       enddo
       return
       end
-      
diff --git a/source/unres/src_MIN/refsys.f b/source/unres/src_MIN/refsys.f
new file mode 100644 (file)
index 0000000..b57c201
--- /dev/null
@@ -0,0 +1,60 @@
+      subroutine refsys(i2,i3,i4,e1,e2,e3,fail)
+      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 (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.
+      logical fail
+      double precision e1(3),e2(3),e3(3)
+      double precision u(3),z(3)
+      include 'COMMON.IOUNITS'
+      include "COMMON.CHAIN"
+      data coinc /1.0d-13/,align /1.0d-13/
+      fail=.false.
+      s1=0.0d0
+      s2=0.0d0
+      do 1 i=1,3
+      zi=c(i,i2)-c(i,i3)
+      ui=c(i,i4)-c(i,i3)
+      s1=s1+zi*zi
+      s2=s2+ui*ui
+      z(i)=zi
+    1 u(i)=ui
+      s1=sqrt(s1)
+      s2=sqrt(s2)
+      if (s1.gt.coinc) goto 2
+      write (iout,1000) i2,i3,i1
+      fail=.true.
+      return
+    2 if (s2.gt.coinc) goto 4
+      write(iout,1000) i3,i4,i1
+      fail=.true.
+      return
+    4 s1=1.0/s1
+      s2=1.0/s2
+      v1=z(2)*u(3)-z(3)*u(2)
+      v2=z(3)*u(1)-z(1)*u(3)
+      v3=z(1)*u(2)-z(2)*u(1)
+      anorm=sqrt(v1*v1+v2*v2+v3*v3)
+      if (anorm.gt.align) goto 6
+      write (iout,1010) i2,i3,i4,i1
+      fail=.true.
+      return
+    6 anorm=1.0/anorm
+      e3(1)=v1*anorm
+      e3(2)=v2*anorm
+      e3(3)=v3*anorm
+      e1(1)=z(1)*s1
+      e1(2)=z(2)*s1
+      e1(3)=z(3)*s1
+      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)
+ 1000 format (/1x,' * * * error - atoms',i4,' and',i4,' coincide.')
+ 1010 format (/1x,' * * * error - atoms',2(i4,2h, ),i4,' form a linear')
+      return
+      end
index 459f395..857aec0 100644 (file)
@@ -212,7 +212,6 @@ c        call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
           do k=1,21
             enetb(k,iii+1,iparm)=energia(k)
           enddo
-#define DEBUG
 #ifdef DEBUG
           write (iout,'(2i5,f10.1,3e15.5)') i,iii,
      &     1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree
@@ -236,7 +235,6 @@ c        call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
         call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
         close(ipdb)
 #endif
-#undef DEBUG
         endif
 
         enddo ! iparm
index 87a3495..431680d 100644 (file)
@@ -96,6 +96,7 @@ C Convert sequence to numeric code
       write(iout,*) 'NNT=',NNT,' NCT=',NCT
 c Read distance restraints
       if (constr_dist.gt.0) then
+        if (refstr) call read_ref_structure(*11)
         call read_dist_constr
         call hpb_partition
       endif
@@ -119,6 +120,7 @@ c Read distance restraints
       endif
       write (iout,'(a)')
       return
+   11 stop "Error reading reference structure"
       end
 c-----------------------------------------------------------------------------
       logical function seq_comp(itypea,itypeb,length)
index eb7e462..c13b16e 100644 (file)
@@ -92,6 +92,8 @@
       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
       write (iout,*) "with_dihed_constr ",with_dihed_constr,
      & " CONSTR_DIST",constr_dist
+      refstr = index(controlcard,'REFSTR').gt.0
+      pdbref = index(controlcard,'PDBREF').gt.0
       call flush(iout)
       return
       end