Merge branch 'devel' into UCGM

Conflicts:
	source/unres/energy.f90

diff --git a/PARAM/tube_carbonano2.parm b/PARAM/tube_carbonano2.parm
new file mode 100644 (file)
index 0000000..ea98cb1
--- /dev/null
+++ b/PARAM/tube_carbonano2.parm
@@ -0,0 +1,25 @@
+1.1970470 5.3667307 0 0 0 0 3.0000000
+1.5539975 5.6438808 0 0 0 0 3.0000000
+1.6679316 5.6689787 0 0 0 0 3.0000000
+1.6606077 5.9381499 0 0 0 0 3.0000000
+1.7428987 5.8625088 0 0 0 0 3.0000000
+1.7310307 5.9950466 0 0 0 0 3.0000000
+1.6322831 5.8318806 0 0 0 0 3.0000000
+1.5348705 5.4955850 0 0 0 0 3.0000000
+1.3603992 5.3937664 0 0 0 0 3.0000000
+1.3228511 5.4371481 0 0 0 0 3.0000000
+1.1970470 5.3667307 0 0 0 0 3.0000000
+1.0325602 5.5439558 0 0 0 0 3.0000000
+0.98513186 5.3780737 0 0 0 0 3.0000000
+0.97556829 5.3995867 0 0 0 0 3.0000000
+0.90197319 5.4184709 0 0 0 0 3.0000000
+0.77024281 5.4679136 0 0 0 0 3.0000000
+0.75456488 5.4686551 0 0 0 0 3.0000000
+1.1983876 5.3466215 0 0 0 0 3.0000000
+0.96779823 5.2968884 0 0 0 0 3.0000000
+0.92065424 5.3752089 0 0 0 0 3.0000000
+1.1218165 5.6721835 0 0 0 0 3.0000000
+1.6679316 5.7029562 0 0 0 0 3.0000000
+1.6606077 5.9355397 0 0 0 0 3.0000000
+1.3228511 5.4343948 0 0 0 0 3.0000000
+1.3228511 5.4343948 0 0 0 0 3.0000000

diff --git a/source/unres/CMakeLists.txt b/source/unres/CMakeLists.txt
index 2a48eda..fd62a89 100644 (file)
--- a/source/unres/CMakeLists.txt
+++ b/source/unres/CMakeLists.txt
@@ -72,6 +72,7 @@ if (Fortran_COMPILER_NAME STREQUAL "ifort")
   set (CMAKE_Fortran_FLAGS_RELEASE " ")
   set (CMAKE_Fortran_FLAGS_DEBUG   "-O0 -g ")
   set(FFLAGS0 "-fpp -c -O3 -ip " ) 
+#  set(FFLAGS0 "-CB -g -ip -fpp" ) 
   set(FFLAGS1 "-fpp -c -O " ) 
   set(FFLAGS2 "-fpp -c -g -CA -CB ")
   set(FFLAGS3 "-fpp -c -g -O0 " )

diff --git a/source/unres/control.F90 b/source/unres/control.F90
index 9b5f441..38ac803 100644 (file)
--- a/source/unres/control.F90
+++ b/source/unres/control.F90
@@ -588,24 +588,24 @@
        iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
             ii=nint_gr(i)+1
             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,&
-       iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
+       iatsc_s,iatsc_e,jj+1,nct_molec(1),nint_gr(i),istart(i,ii),iend(i,ii),*12)
 #else
             nint_gr(i)=2
             istart(i,1)=i+1
             iend(i,1)=jj-1
             istart(i,2)=jj+1
-            iend(i,2)=nct
+            iend(i,2)=nct_molec(1)
 #endif
           endif
         else
 #ifdef MPI
           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,&
-          iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
+          iatsc_s,iatsc_e,i+1,nct_molec(1),nint_gr(i),istart(i,1),iend(i,1),*12)
 #else
           nint_gr(i)=1
           istart(i,1)=i+1
-          iend(i,1)=nct
-          ind_scint=ind_scint+nct-i
+          iend(i,1)=nct_molec(1)
+          ind_scint=ind_scint+nct_molec(1)-i
 #endif
         endif
 #ifdef MPI
@@ -648,16 +648,16 @@
       iatel_e=0
       ind_eleint=0
       ind_eleint_old=0
-      if (itype(nres_molec(1),1).eq.ntyp_molec(1)) then
+      if (itype(nres_molec(1),1).eq.ntyp1_molec(1)) then
       nct_molec(1)=nres_molec(1)-1
       else
       nct_molec(1)=nres_molec(1)
       endif
-
+!       print *,"nct",nct,nct_molec(1),itype(nres_molec(1),1),ntyp_molec(1)
       do i=nnt,nct_molec(1)-3
         ijunk=0
         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,&
-          iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
+          iatel_s,iatel_e,i+ispp,nct_molec(1)-1,ijunk,ielstart(i),ielend(i),*13)
       enddo ! i 
    13 continue
       if (iatel_s.eq.0) iatel_s=1
@@ -674,7 +674,7 @@
         ijunk=0
         call int_partition(ind_eleint_vdw,my_ele_inds_vdw,&
           my_ele_inde_vdw,i,&
-          iatel_s_vdw,iatel_e_vdw,i+2,nct-1,ijunk,ielstart_vdw(i),&
+          iatel_s_vdw,iatel_e_vdw,i+2,nct_molec(1)-1,ijunk,ielstart_vdw(i),&
           ielend_vdw(i),*15)
 !        write (iout,*) i," ielstart_vdw",ielstart_vdw(i),
 !     &   " ielend_vdw",ielend_vdw(i)
@@ -721,7 +721,7 @@
         if (i.lt.nnt+iscp) then
 !d        write (iout,*) 'i.le.nnt+iscp'
           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,&
-            iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),&
+            iatscp_s,iatscp_e,i+iscp,nct_molec(1),nscp_gr(i),iscpstart(i,1),&
             iscpend(i,1),*14)
         else if (i.gt.nct-iscp) then
 !d        write (iout,*) 'i.gt.nct-iscp'
@@ -734,7 +734,7 @@
            iscpend(i,1),*14)
           ii=nscp_gr(i)+1
           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,&
-            iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),&
+            iatscp_s,iatscp_e,i+iscp,nct_molec(1),nscp_gr(i),iscpstart(i,ii),&
             iscpend(i,ii),*14)
         endif
       enddo ! i
@@ -796,10 +796,11 @@
       iphid_end=iturn4_end+2
       iturn4_start=iturn4_start-1
       iturn4_end=iturn4_end-1
-      print *,"TUTUTU",nres_molec(1),nres
-      call int_bounds(nres-2,ibond_start,ibond_end) 
+!      print *,"TUTUTU",nres_molec(1),nres
+      call int_bounds(nres_molec(1)-2,ibond_start,ibond_end) 
       ibond_start=ibond_start+1
       ibond_end=ibond_end+1
+      print *,ibond_start,ibond_end
       call int_bounds(nct_molec(1)-nnt,ibondp_start,ibondp_end) 
       ibondp_start=ibondp_start+nnt
       ibondp_end=ibondp_end+nnt

diff --git a/source/unres/energy.f90 b/source/unres/energy.f90
index c21d21d..3975993 100644 (file)
--- a/source/unres/energy.f90
+++ b/source/unres/energy.f90
@@ -45,6 +45,7 @@
       real(kind=8),dimension(:,:,:),allocatable :: gacont      !(3,maxconts,maxres)
       integer,dimension(:),allocatable :: ishield_list
       integer,dimension(:,:),allocatable ::  shield_list
+      real(kind=8),dimension(:),allocatable :: enetube,enecavtube
 !                
 ! 12/26/95 - H-bonding contacts
 !      common /contacts_hb/ 
@@ -421,6 +422,7 @@
 ! Calculate the bond-stretching energy
 !
       call ebond(estr)
+       print *,"EBOND",estr
 !       write(iout,*) "in etotal afer ebond",ipot
 
 ! 
@@ -5143,11 +5145,13 @@
 !        endif
       enddo
       estr=0.5d0*AKP*estr+estr1
+      print *,"estr_bb",estr,AKP
 !
 ! 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
 !
       do i=ibond_start,ibond_end
         iti=iabs(itype(i,1))
+        if (iti.eq.0) print *,"WARNING WRONG SETTTING",i
         if (iti.ne.10 .and. iti.ne.ntyp1) then
           nbi=nbondterm(iti)
           if (nbi.eq.1) then
@@ -5156,6 +5160,7 @@
             "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,&
             AKSC(1,iti),AKSC(1,iti)*diff*diff
             estr=estr+0.5d0*AKSC(1,iti)*diff*diff
+            print *,"estr_sc",estr
             do j=1,3
               gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
             enddo
@@ -5184,7 +5189,9 @@
               usumsqder=usumsqder+ud(j)*uprod2   
             enddo
             estr=estr+uprod/usum
-            if (energy_dec) write (iout,*) &
+            print *,"estr_sc",estr,i
+
+             if (energy_dec) write (iout,*) &
             "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,&
             AKSC(1,iti),AKSC(1,iti)*diff*diff
             do j=1,3
@@ -18330,7 +18337,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !C and r0 is the excluded size of nanotube (can be set to 0 if we want just a 
 !C simple Kihara potential
       subroutine calctube(Etube)
-      real(kind=8) :: vectube(3),enetube(nres*2)
+      real(kind=8),dimension(3) :: vectube
       real(kind=8) :: Etube,xtemp,xminact,yminact,& 
        ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi, &
        sc_aa_tube,sc_bb_tube
@@ -18491,7 +18498,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !C and r0 is the excluded size of nanotube (can be set to 0 if we want just a 
 !C simple Kihara potential
       subroutine calctube2(Etube)
-      real(kind=8) :: vectube(3),enetube(nres*2)
+            real(kind=8),dimension(3) :: vectube
       real(kind=8) :: Etube,xtemp,xminact,yminact,&
        ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi,fracinbuf,&
        sstube,ssgradtube,sc_aa_tube,sc_bb_tube
@@ -18728,8 +18735,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         end subroutine calctube2
 !=====================================================================================================================================
       subroutine calcnano(Etube)
-      real(kind=8) :: vectube(3),enetube(nres*2), &
-      enecavtube(nres*2)
+      real(kind=8),dimension(3) :: vectube
+      
       real(kind=8) :: Etube,xtemp,xminact,yminact,&
        ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,denominator,faccav,&
        sc_aa_tube,sc_bb_tube,zmin,ztemp,zminact
@@ -18829,7 +18836,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !C         fac=fac+faccav
 !C 667     continue
          endif
-
+          if (energy_dec) write(iout,*),i,rdiff,enetube(i),enecavtube(i)
         do j=1,3
         gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0
         gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0
@@ -18932,6 +18939,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac
           gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
          enddo
+          if (energy_dec) write(iout,*),i,rdiff,enetube(i+nres),enecavtube(i+nres)
         enddo
 
 
@@ -19555,6 +19563,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
       allocate(shield_list(50,nres))
       allocate(dyn_ss_mask(nres))
       allocate(fac_shield(nres))
+      allocate(enetube(nres*2))
+      allocate(enecavtube(nres*2))
+
 !(maxres)
       dyn_ss_mask(:)=.false.
 !----------------------

diff --git a/source/unres/geometry.f90 b/source/unres/geometry.f90
index 7b6c961..8542f52 100644 (file)
--- a/source/unres/geometry.f90
+++ b/source/unres/geometry.f90
@@ -3039,7 +3039,7 @@
       integer :: i,j,ires,nscat
       real(kind=8),dimension(3,20) :: sccor
       real(kind=8) :: sccmj
-        print *,"I am in sccenter",ires,nscat
+!        print *,"I am in sccenter",ires,nscat
       do j=1,3
         sccmj=0.0D0
         do i=1,nscat

diff --git a/source/unres/io.f90 b/source/unres/io.f90
index 24549e8..ae6e636 100644 (file)
--- a/source/unres/io.f90
+++ b/source/unres/io.f90
@@ -720,7 +720,7 @@
 !      integer :: rescode
 !      double precision x(maxvar)
       character(len=256) :: pdbfile
-      character(len=320) :: weightcard
+      character(len=480) :: weightcard
       character(len=80) :: weightcard_t!,ucase
 !      integer,dimension(:),allocatable :: itype_pdb   !(maxres)
 !      common /pizda/ itype_pdb

diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90
index cf23282..5153cc9 100644 (file)
--- a/source/unres/io_config.f90
+++ b/source/unres/io_config.f90
@@ -2466,7 +2466,7 @@
               ires=ires-ishift+ishift1
               ires_old=ires
             endif 
-            print *,'atom',ires,atom
+!            print *,'atom',ires,atom
             if (res.eq.'ACE' .or. res.eq.'NHE') then
               itype(ires,1)=10
             else
@@ -2524,7 +2524,7 @@
 !            enddo
              counter=0
            endif
-            print *, "ATOM",atom(1:3)
+!            print *, "ATOM",atom(1:3)
           else if (atom(1:3).eq."C5'") then
              read (card(19:19),'(a1)') sugar
              isugar=sugarcode(sugar,ires)
@@ -2654,7 +2654,7 @@
 !      nres=ires
       nsup=nres
       nstart_sup=1
-      print *,"molecule",molecule
+!      print *,"molecule",molecule
       if (itype(nres,1).ne.10) then
         nres=nres+1
         itype(nres,molecule)=ntyp1_molec(molecule)
@@ -2678,7 +2678,7 @@
         enddo
         endif
       endif
-     print *,'I have',nres, nres_molec(:)
+!     print *,'I have',nres, nres_molec(:)
 
 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
 #ifdef WHAM_RUN
@@ -2825,6 +2825,29 @@
          itype(i,k)=itype_temporary(i,k)
         enddo
        enddo
+      if (itype(1,1).eq.ntyp1) then
+        nsup=nsup-1
+        nstart_sup=2
+        if (unres_pdb) then
+! 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)-1.9d0*e2(j)
+          enddo
+        else
+        do j=1,3
+          dcj=(c(j,4)-c(j,3))/2.0
+          c(j,1)=c(j,2)-dcj
+          c(j,nres+1)=c(j,1)
+        enddo
+        endif
+      endif
+
       if (lprn) then
       write (iout,'(/a)') &
         "Cartesian coordinates of the reference structure after sorting"