pierwsza proba merga devela do adasko
authorAdam Sieradzan <adasko@piasek2.(none)>
Tue, 5 Nov 2013 13:37:15 +0000 (14:37 +0100)
committerAdam Sieradzan <adasko@piasek2.(none)>
Tue, 5 Nov 2013 13:37:15 +0000 (14:37 +0100)
Merge branch 'devel' into adasko

Conflicts:
.gitignore
PARAM/pot_theta_G631_DIL.parm
PARAM/pot_theta_G631_DIL_ext.parm
bin/unres/MD-M/unres_Tc_procor_oldparm_em64-D-symetr.exe
bin/unres_clustMD_MPI-oldparm
bin/wham/wham_multparm-ham_rep-oldparm
source/cluster/wham/src-M/CMakeLists.txt
source/cluster/wham/src-M/COMMON.SCCOR
source/cluster/wham/src-M/cartprint.o
source/cluster/wham/src-M/chainbuild.o
source/cluster/wham/src-M/contact.o
source/cluster/wham/src-M/convert.o
source/cluster/wham/src-M/energy_p_new.o
source/cluster/wham/src-M/geomout.o
source/cluster/wham/src-M/hc.o
source/cluster/wham/src-M/initialize_p.o
source/cluster/wham/src-M/int_from_cart1.o
source/cluster/wham/src-M/main_clust.o
source/cluster/wham/src-M/parmread.o
source/cluster/wham/src-M/permut.o
source/cluster/wham/src-M/proc_proc.o
source/cluster/wham/src-M/read_coords.o
source/cluster/wham/src-M/read_ref_str.o
source/cluster/wham/src-M/readpdb.o
source/cluster/wham/src-M/readrtns.o
source/cluster/wham/src-M/rescode.o
source/cluster/wham/src-M/setup_var.o
source/cluster/wham/src-M/wrtclust.o
source/cluster/wham/src-M/xdrf/libxdrf.a
source/cluster/wham/src/energy_p_new.F
source/unres/src_MD-M/COMMON.SCCOR
source/unres/src_MD-M/COMMON.TORSION
source/unres/src_MD-M/Makefile
source/unres/src_MD-M/cinfo.f
source/unres/src_MD-M/stochfric.F
source/unres/src_MD/COMMON.LOCAL
source/unres/src_MD/COMMON.SCCOR
source/unres/src_MD/Makefile
source/unres/src_MD/cinfo.f
source/unres/src_MD/parmread.F
source/wham/src-M/CMakeLists.txt
source/wham/src-M/cinfo.f
source/wham/src-M/include_unres/COMMON.SCCOR
source/wham/src/Makefile
source/wham/src/cinfo.f
source/wham/src/enecalc1.F
source/wham/src/energy_p_new.F
source/wham/src/include_unres/COMMON.SCCOR
source/wham/src/int_from_cart.f
source/wham/src/parmread.F

32 files changed:
1  2 
.gitignore
source/cluster/wham/src-M/energy_p_new.F
source/cluster/wham/src/energy_p_new.F
source/cluster/wham/src/initialize_p.F
source/cluster/wham/src/parmread.F
source/cluster/wham/src/read_coords.F
source/cluster/wham/src/readrtns.F
source/unres/src_MD-M/CMakeLists.txt
source/unres/src_MD-M/COMMON.CHAIN
source/unres/src_MD-M/MD_A-MTS.F
source/unres/src_MD-M/Makefile
source/unres/src_MD-M/energy_p_new_barrier.F
source/unres/src_MD-M/initialize_p.F
source/unres/src_MD-M/readpdb.F
source/unres/src_MD-M/stochfric.F
source/unres/src_MD/CMakeLists.txt
source/unres/src_MD/COMMON.SCCOR
source/unres/src_MD/energy_p_new_barrier.F
source/unres/src_MD/parmread.F
source/unres/src_MD/readpdb.F
source/unres/src_MD/readrtns.F
source/unres/src_MD/sc_move.F
source/unres/src_MD/stochfric.F
source/wham/src-M/CMakeLists.txt
source/wham/src-M/energy_p_new.F
source/wham/src/COMMON.ALLPARM
source/wham/src/enecalc1.F
source/wham/src/energy_p_new.F
source/wham/src/initialize_p.F
source/wham/src/molread_zs.F
source/wham/src/parmread.F
source/wham/src/store_parm.F

diff --cc .gitignore
@@@ -10,14 -11,17 +11,19 @@@ cinfo.
  
  # ignore build dir
  build/
 +build2/
  
+ # latex files in documentation 
+ doc/*/*.aux
+ doc/*/*.log
  # ignored dirs form adasko
  gradcheck/
  mapcheck/
  run/
  sympcheck/
 -compinfo
 -DIL/
  bin/unres/MD/unres_ifort_MPICH_GAB_czyt.exe
 +bin/unres/MD-M/unres_Tc_procor_newparm_em64-D-symetr.exe
- DIL
++DIL/
 +compinfo
@@@ -748,8 -750,14 +750,14 @@@ C Calculate SC interaction energy
  C
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
+             IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+               call dyn_ssbond_ene(i,j,evdwij)
+               evdw=evdw+evdwij
+ c              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+ c     &                        'evdw',i,j,evdwij,' ss'
+             ELSE
              ind=ind+1
 -            itypj=itype(j)
 +            itypj=iabs(itype(j))
              dscj_inv=vbld_inv(j+nres)
              sig0ij=sigma(itypi,itypj)
              chi1=chi(itypi,itypj)
@@@ -879,8 -889,15 +889,15 @@@ C Calculate SC interaction energy
  C
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
+ C in case of diagnostics    write (iout,*) "TU SZUKAJ",i,j,dyn_ss_mask(i),dyn_ss_mask(j)
+             IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+               call dyn_ssbond_ene(i,j,evdwij)
+               evdw=evdw+evdwij
+ c              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+ c     &                        'evdw',i,j,evdwij,' ss'
+             ELSE
              ind=ind+1
 -            itypj=itype(j)
 +            itypj=iabs(itype(j))
              dscj_inv=vbld_inv(j+nres)
              sig0ij=sigma(itypi,itypj)
              r0ij=r0(itypi,itypj)
@@@ -2822,8 -2840,9 +2840,10 @@@ c        write (iout,*) "i",i," ii",ii,
  c     &    dhpb(i),dhpb1(i),forcon(i)
  C 24/11/03 AL: SS bridges handled separately because of introducing a specific
  C    distance and angle dependent SS bond potential.
+         if (.not.dyn_ss .and. i.le.nss) then
+ C 15/02/13 CC dynamic SSbond - additional check
 -        if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
 +        if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
 +     &  iabs( itype(jjj)).eq.1) then
            call ssbond_ene(iii,jjj,eij)
            ehpb=ehpb+2*eij
  cd          write (iout,*) "eij",eij
Simple merge
Simple merge
Simple merge
Simple merge
Simple merge
Simple merge
Simple merge
Simple merge
@@@ -138,11 -140,18 +140,18 @@@ C Calculate the CM of the last side cha
        nstart_sup=1
        if (itype(nres).ne.10) then
          nres=nres+1
 -        itype(nres)=21
 +        itype(nres)=ntyp1
          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)
        include 'COMMON.INTERACT'
        include 'COMMON.IOUNITS'
        include 'COMMON.NAMES'
-       double precision radius(maxres2),gamvec(maxres6)
+       double precision radius(maxres2),gamvec(maxres2)
        parameter (twosix=1.122462048309372981d0)
 -      logical lprn /.true./
 +      logical lprn /.false./
  c
  c     determine new friction coefficients every few SD steps
  c
Simple merge
Simple merge
@@@ -5806,15 -5786,8 +5797,10 @@@ c     lprn=.true
        etors=0.0D0
        do i=iphi_start,iphi_end
        etors_ii=0.0D0
 +c        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
 +c     &      .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
          itori=itortyp(itype(i-2))
          itori1=itortyp(itype(i-1))
-         if (iabs(itype(i)).eq.20) then
-         iblock=2
-         else
-         iblock=1
-         endif
          phii=phi(i)
          gloci=0.0D0
  C Regular cosine and sine terms
@@@ -673,16 -490,12 +605,20 @@@ c      write (iout,*) 'ntortyp',ntorty
  C
  C 6/23/01 Read parameters for double torsionals
  C
-       do iblock=1,2
-       do i=0,ntortyp-1
-         do j=-ntortyp+1,ntortyp-1
-           do k=-ntortyp+1,ntortyp-1
+       do i=1,ntortyp
+         do j=1,ntortyp
+           do k=1,ntortyp
              read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
++<<<<<<< HEAD
 +c              write (iout,*) "OK onelett",
 +c     &         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
++>>>>>>> devel
                write (iout,*) "Error in double torsional parameter file",
       &         i,j,k,t1,t2,t3
  #ifdef MPI
        if (lprint) then
        write (iout,*) 
        write (iout,*) 'Constants for double torsionals'
++<<<<<<< HEAD
 +      do iblock=1,2
 +      do i=0,ntortyp-1
 +        do j=-ntortyp+1,ntortyp-1
 +          do k=-ntortyp+1,ntortyp-1
++=======
+       do i=1,ntortyp
+         do j=1,ntortyp 
+           do k=1,ntortyp
++>>>>>>> devel
              write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
-      &        ' nsingle',ntermd_1(i,j,k,iblock),
-      &        ' ndouble',ntermd_2(i,j,k,iblock)
+      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
              write (iout,*)
              write (iout,*) 'Single angles:'
-             do l=1,ntermd_1(i,j,k,iblock)
-               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
-      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
-      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
-      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
+             do l=1,ntermd_1(i,j,k)
+               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
+      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
+      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
              enddo
              write (iout,*)
              write (iout,*) 'Pairs of angles:'
@@@ -771,71 -557,27 +687,90 @@@ C Read of Side-chain backbone correlati
  C Modified 11 May 2012 by Adasko
  CCC
  C
++<<<<<<< HEAD
 +      read (isccor,*,end=119,err=119) nsccortyp
 +#ifdef SCCORPDB
 +      read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 +      do i=-ntyp,-1
 +        isccortyp(i)=-isccortyp(-i)
 +      enddo
 +      iscprol=isccortyp(20)
++=======
+       read (isccor,*,end=1113,err=1113) nsccortyp
+       read (isccor,*,end=1113,err=1113) (isccortyp(i),i=1,ntyp)
++>>>>>>> devel
  c      write (iout,*) 'ntortyp',ntortyp
        maxinter=3
  cc maxinter is maximum interaction sites
        do l=1,maxinter    
        do i=1,nsccortyp
        do j=1,nsccortyp
++<<<<<<< HEAD
 +        read (isccor,*,end=119,err=119) nterm_sccor(i,j),nlor_sccor(i,j)
++=======
+         read (isccor,*,end=1113,err=1113) nterm_sccor(i,j),
+      &             nlor_sccor(i,j)
++>>>>>>> devel
            v0ijsccor=0.0d0
 +          v0ijsccor1=0.0d0
 +          v0ijsccor2=0.0d0
 +          v0ijsccor3=0.0d0
            si=-1.0d0
 -  
 +          nterm_sccor(-i,j)=nterm_sccor(i,j)
 +          nterm_sccor(-i,-j)=nterm_sccor(i,j)
 +          nterm_sccor(i,-j)=nterm_sccor(i,j)  
          do k=1,nterm_sccor(i,j)
++<<<<<<< HEAD
 +          read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 +     &    ,v2sccor(k,l,i,j)
 +            if (j.eq.iscprol) then
 +              if (i.eq.isccortyp(10)) then
 +              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 +              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 +              else
 +             v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 +     &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 +             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 +     &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 +             v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 +             v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 +             v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 +             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)          
 +             endif
 +            else
 +              if (i.eq.isccortyp(10)) then
 +              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 +              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 +              else
 +                if (j.eq.isccortyp(10)) then
 +              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
 +              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
 +                else
 +             v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
 +             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 +             v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 +             v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 +             v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 +             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 +            endif
 +             endif
 +             endif 
++=======
+           read (isccor,*,end=1113,err=1113) kk,v1sccor(k,l,i,j)
+      &    ,v2sccor(k,l,i,j) 
++>>>>>>> devel
              v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 +            v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 +            v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 +            v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
              si=-si
            enddo
          do k=1,nlor_sccor(i,j)
++<<<<<<< HEAD
 +            read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
++=======
+             read (isccor,*,end=1113,err=1113) kk,vlor1sccor(k,i,j),
++>>>>>>> devel
       &        vlor2sccor(k,i,j),vlor3sccor(k,i,j) 
              v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
       &(1+vlor3sccor(k,i,j)**2)
  c        b1(1,i)=0.0d0
  c        b1(2,i)=0.0d0
          B1tilde(1,i) = b(3)
++<<<<<<< HEAD
 +        B1tilde(2,i) =-b(5)
 +        B1tilde(1,-i) =-b(3)
 +        B1tilde(2,-i) =b(5)
++=======
+         B1tilde(2,i) =-b(5) 
++>>>>>>> devel
  c        b1tilde(1,i)=0.0d0
  c        b1tilde(2,i)=0.0d0
          B2(1,i)  = b(2)
@@@ -102,11 -104,18 +104,18 @@@ C Calculate the CM of the last side cha
        nstart_sup=1
        if (itype(nres).ne.10) then
          nres=nres+1
 -        itype(nres)=21
 +        itype(nres)=ntyp1
          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)
Simple merge
Simple merge
Simple merge
@@@ -112,20 -112,24 +112,36 @@@ endif(UNRES_WITH_MPI
  set_property(SOURCE ${UNRES_WHAM_M_SRC0} PROPERTY COMPILE_FLAGS ${FFLAGS0} )
  
  #=========================================
- # WHAM preprocesor flags
+ #  Settings for GAB force field
  #=========================================
+ if(UNRES_MD_FF STREQUAL "GAB" )
+   # set preprocesor flags   
+   set(CPPFLAGS "PROCOR  -DSPLITELE -DCRYST_BOND -DCRYST_THETA -DCRYST_SC  -DSCCORPDB" )
  
++<<<<<<< HEAD
 +if(UNRES_MD_FF STREQUAL "GAB" )
 +  # set preprocesor flags   
 +  set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC  -DSCCORPDB" )
 +
++=======
++>>>>>>> devel
  #=========================================
  #  Settings for E0LL2Y force field
  #=========================================
  elseif(UNRES_MD_FF STREQUAL "E0LL2Y")
    # set preprocesor flags   
++<<<<<<< HEAD
 +  set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DSCCORPDB" )
 +endif(UNRES_MD_FF STREQUAL "GAB")
++=======
+   set(CPPFLAGS "PROCOR  -DSPLITELE -DSCCORPDB" )
+ endif(UNRES_MD_FF STREQUAL "GAB")
+ #=========================================
+ # Additional flags
+ #=========================================
+ set(CPPFLAGS "${CPPFLAGS} -DUNRES -DISNAN") 
++>>>>>>> devel
  
  #=========================================
  # System specific flags
Simple merge
Simple merge
@@@ -194,12 -195,12 +195,16 @@@ c        call pdbout(ii+1,beta_h(ib,ipa
       &         " the value read in: ",energia(0),eini," point",
       &         iii+1,indstart(me1)+iii," T",
       &         1.0d0/(1.987D-3*beta_h(ib,ipar))
-               errmsg_count=errmsg_count+1
+              call enerprint(energia(0),fT)
+              call pdbout(iii+1,beta_h(ib,ipar),
+      &                   eini,energia(0),0.0d0,rmsdev)
+              write (iout,*)
  
+               errmsg_count=errmsg_count+1
 +              call pdbout(indstart(me1)+iii,
 +     & 1.0d0/(1.987D-3*beta_h(ib,ipar)),
 +     &energia(0),eini,0.0d0,0.0d0)
 +              call enerprint(energia(0),fT)
                if (errmsg_count.gt.maxerrmsg_count) 
       &          write (iout,*) "Too many warning messages"
                if (einicheck.gt.1) then
@@@ -804,8 -802,23 +806,23 @@@ C Calculate SC interaction energy
  C
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
+ C in case of diagnostics    write (iout,*) "TU SZUKAJ",i,j,dyn_ss_mask(i),dyn_ss_mask(j)
+ C /06/28/2013 Adasko: In case of dyn_ss - dynamic disulfide bond
+ C formation no electrostatic interactions should be calculated. If it
+ C would be allowed NaN would appear
+             IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+ C /06/28/2013 Adasko: dyn_ss_mask is logical statement wheather this Cys
+ C residue can or cannot form disulfide bond. There is still bug allowing
+ C Cys...Cys...Cys bond formation
+               call dyn_ssbond_ene(i,j,evdwij)
+ C /06/28/2013 Adasko: dyn_ssbond_ene is dynamic SS bond foration energy
+ C function in ssMD.F
+               evdw=evdw+evdwij
+ c              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+ c     &                        'evdw',i,j,evdwij,' ss'
+             ELSE
              ind=ind+1
 -            itypj=itype(j)
 +            itypj=iabs(itype(j))
              dscj_inv=vbld_inv(j+nres)
              sig0ij=sigma(itypi,itypj)
              chi1=chi(itypi,itypj)
@@@ -2901,10 -2916,12 +2920,13 @@@ c        write (iout,*) "i",i," ii",ii,
  c     &    dhpb(i),dhpb1(i),forcon(i)
  C 24/11/03 AL: SS bridges handled separately because of introducing a specific
  C    distance and angle dependent SS bond potential.
-         if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
+         if (.not.dyn_ss .and. i.le.nss) then
+ C 15/02/13 CC dynamic SSbond - additional check
 -        if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
++         if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
 +     & iabs(itype(jjj)).eq.1) then
            call ssbond_ene(iii,jjj,eij)
            ehpb=ehpb+2*eij
+          endif
  cd          write (iout,*) "eij",eij
          else if (ii.gt.nres .and. jj.gt.nres) then
  c Restraints from contact prediction
Simple merge
Simple merge
@@@ -640,13 -551,12 +662,13 @@@ c      write (iout,*) 'ntortyp',ntorty
  C
  C 6/23/01 Read parameters for double torsionals
  C
 -      do i=1,ntortyp
 -        do j=1,ntortyp
 -          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)') 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
                 stop "Error in double torsional parameter file"
Simple merge