character(len=1) :: t1,t2,t3
character(len=1) :: onelett(4) = (/"G","A","P","D"/)
character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
- logical :: lprint,LaTeX
+ logical :: lprint,LaTeX,SPLIT_FOURIERTOR
real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
real(kind=8),dimension(13) :: buse
character(len=3) :: lancuch !,ucase
!el local variables
- integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
- integer :: maxinter,junk,kk,ii,ncatprotparm
+ integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj
+ integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp
real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
sigmasctube
- integer :: ichir1,ichir2
+ integer :: ichir1,ichir2,ijunk
+ character*3 string
+
! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
!el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
!el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
! Read the parameters of Utheta determined from ab initio surfaces
! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
!
+ IF (tor_mode.eq.0) THEN
read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
ntheterm3,nsingle,ndouble
nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
enddo
call flush(iout)
endif
+ ELSE
+!C here will be the apropriate recalibrating for D-aminoacid
+ read (ithep,*,end=121,err=121) nthetyp
+ allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
+ allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
+ do i=-nthetyp+1,nthetyp-1
+ read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
+ do j=0,nbend_kcc_Tb(i)
+ read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
+ enddo
+ enddo
+ if (lprint) then
+ write (iout,'(a)') &
+ "Parameters of the valence-only potentials"
+ do i=-nthetyp+1,nthetyp-1
+ write (iout,'(2a)') "Type ",toronelet(i)
+ do k=0,nbend_kcc_Tb(i)
+ write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
+ enddo
+ enddo
+ endif
+ ENDIF ! TOR_MODE
+
write (2,*) "Start reading THETA_PDB",ithep_pdb
do i=1,ntyp
! write (2,*) 'i=',i
#endif
close(irotam)
+
+!C
+!C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
+!C interaction energy of the Gly, Ala, and Pro prototypes.
+!C
+ read (ifourier,*) nloctyp
+ SPLIT_FOURIERTOR = nloctyp.lt.0
+ nloctyp = iabs(nloctyp)
+!C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
+!C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
+!C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
+!C allocate(ctilde(2,2,nres))
+!C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
+!C allocate(gtb1(2,nres))
+!C allocate(gtb2(2,nres))
+!C allocate(cc(2,2,nres))
+!C allocate(dd(2,2,nres))
+!C allocate(ee(2,2,nres))
+!C allocate(gtcc(2,2,nres))
+!C allocate(gtdd(2,2,nres))
+!C allocate(gtee(2,2,nres))
+
+#ifdef NEWCORR
+ allocate(itype2loc(-ntyp1:ntyp1))
+ allocate(iloctyp(-nloctyp:nloctyp))
+ allocate(bnew1(3,2,-nloctyp:nloctyp))
+ allocate(bnew2(3,2,-nloctyp:nloctyp))
+ allocate(ccnew(3,2,-nloctyp:nloctyp))
+ allocate(ddnew(3,2,-nloctyp:nloctyp))
+ allocate(e0new(3,-nloctyp:nloctyp))
+ allocate(eenew(2,2,2,-nloctyp:nloctyp))
+ allocate(bnew1tor(3,2,-nloctyp:nloctyp))
+ allocate(bnew2tor(3,2,-nloctyp:nloctyp))
+ allocate(ccnewtor(3,2,-nloctyp:nloctyp))
+ allocate(ddnewtor(3,2,-nloctyp:nloctyp))
+ allocate(e0newtor(3,-nloctyp:nloctyp))
+ allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
+
+ read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
+ read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
+ itype2loc(ntyp1)=nloctyp
+ iloctyp(nloctyp)=ntyp1
+ do i=1,ntyp1
+ itype2loc(-i)=-itype2loc(i)
+ enddo
+#else
+ allocate(iloctyp(-nloctyp:nloctyp))
+ allocate(itype2loc(-ntyp1:ntyp1))
+ iloctyp(0)=10
+ iloctyp(1)=9
+ iloctyp(2)=20
+ iloctyp(3)=ntyp1
+#endif
+ do i=1,nloctyp
+ iloctyp(-i)=-iloctyp(i)
+ enddo
+!c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
+!c write (iout,*) "nloctyp",nloctyp,
+!c & " iloctyp",(iloctyp(i),i=0,nloctyp)
+!c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
+!c write (iout,*) "nloctyp",nloctyp,
+!c & " iloctyp",(iloctyp(i),i=0,nloctyp)
+#ifdef NEWCORR
+ do i=0,nloctyp-1
+!c write (iout,*) "NEWCORR",i
+ read (ifourier,*,end=115,err=115)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR BNEW1"
+!c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR BNEW2"
+!c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
+ read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
+ enddo
+!c write (iout,*) "NEWCORR CCNEW"
+!c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
+ read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
+ enddo
+!c write (iout,*) "NEWCORR DDNEW"
+!c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,2
+ do jj=1,2
+ do kk=1,2
+ read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
+ enddo
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR EENEW1"
+!c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
+ do ii=1,3
+ read (ifourier,*,end=115,err=115) e0new(ii,i)
+ enddo
+!c write (iout,*) (e0new(ii,i),ii=1,3)
+ enddo
+!c write (iout,*) "NEWCORR EENEW"
+ do i=0,nloctyp-1
+ do ii=1,3
+ ccnew(ii,1,i)=ccnew(ii,1,i)/2
+ ccnew(ii,2,i)=ccnew(ii,2,i)/2
+ ddnew(ii,1,i)=ddnew(ii,1,i)/2
+ ddnew(ii,2,i)=ddnew(ii,2,i)/2
+ enddo
+ enddo
+ do i=1,nloctyp-1
+ do ii=1,3
+ bnew1(ii,1,-i)= bnew1(ii,1,i)
+ bnew1(ii,2,-i)=-bnew1(ii,2,i)
+ bnew2(ii,1,-i)= bnew2(ii,1,i)
+ bnew2(ii,2,-i)=-bnew2(ii,2,i)
+ enddo
+ do ii=1,3
+!c ccnew(ii,1,i)=ccnew(ii,1,i)/2
+!c ccnew(ii,2,i)=ccnew(ii,2,i)/2
+!c ddnew(ii,1,i)=ddnew(ii,1,i)/2
+!c ddnew(ii,2,i)=ddnew(ii,2,i)/2
+ ccnew(ii,1,-i)=ccnew(ii,1,i)
+ ccnew(ii,2,-i)=-ccnew(ii,2,i)
+ ddnew(ii,1,-i)=ddnew(ii,1,i)
+ ddnew(ii,2,-i)=-ddnew(ii,2,i)
+ enddo
+ e0new(1,-i)= e0new(1,i)
+ e0new(2,-i)=-e0new(2,i)
+ e0new(3,-i)=-e0new(3,i)
+ do kk=1,2
+ eenew(kk,1,1,-i)= eenew(kk,1,1,i)
+ eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
+ eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
+ eenew(kk,2,2,-i)= eenew(kk,2,2,i)
+ enddo
+ enddo
+ if (lprint) then
+ write (iout,'(a)') "Coefficients of the multibody terms"
+ do i=-nloctyp+1,nloctyp-1
+ write (iout,*) "Type: ",onelet(iloctyp(i))
+ write (iout,*) "Coefficients of the expansion of B1"
+ do j=1,2
+ write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of B2"
+ do j=1,2
+ write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of C"
+ write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
+ write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of D"
+ write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
+ write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of E"
+ write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
+ do j=1,2
+ do k=1,2
+ write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
+ enddo
+ enddo
+ enddo
+ endif
+ IF (SPLIT_FOURIERTOR) THEN
+ do i=0,nloctyp-1
+!c write (iout,*) "NEWCORR TOR",i
+ read (ifourier,*,end=115,err=115)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR BNEW1 TOR"
+!c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR BNEW2 TOR"
+!c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
+ read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
+ enddo
+!c write (iout,*) "NEWCORR CCNEW TOR"
+!c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
+ read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
+ enddo
+!c write (iout,*) "NEWCORR DDNEW TOR"
+!c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,2
+ do jj=1,2
+ do kk=1,2
+ read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
+ enddo
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR EENEW1 TOR"
+!c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
+ do ii=1,3
+ read (ifourier,*,end=115,err=115) e0newtor(ii,i)
+ enddo
+!c write (iout,*) (e0newtor(ii,i),ii=1,3)
+ enddo
+!c write (iout,*) "NEWCORR EENEW TOR"
+ do i=0,nloctyp-1
+ do ii=1,3
+ ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
+ ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
+ ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
+ ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
+ enddo
+ enddo
+ do i=1,nloctyp-1
+ do ii=1,3
+ bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
+ bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
+ bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
+ bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
+ enddo
+ do ii=1,3
+ ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
+ ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
+ ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
+ ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
+ enddo
+ e0newtor(1,-i)= e0newtor(1,i)
+ e0newtor(2,-i)=-e0newtor(2,i)
+ e0newtor(3,-i)=-e0newtor(3,i)
+ do kk=1,2
+ eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
+ eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
+ eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
+ eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
+ enddo
+ enddo
+ if (lprint) then
+ write (iout,'(a)') &
+ "Single-body coefficients of the torsional potentials"
+ do i=-nloctyp+1,nloctyp-1
+ write (iout,*) "Type: ",onelet(iloctyp(i))
+ write (iout,*) "Coefficients of the expansion of B1tor"
+ do j=1,2
+ write (iout,'(3hB1(,i1,1h),3f10.5)') &
+ j,(bnew1tor(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of B2tor"
+ do j=1,2
+ write (iout,'(3hB2(,i1,1h),3f10.5)') &
+ j,(bnew2tor(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of Ctor"
+ write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
+ write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of Dtor"
+ write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
+ write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of Etor"
+ write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
+ do j=1,2
+ do k=1,2
+ write (iout,'(1hE,2i1,2f10.5)') &
+ j,k,(eenewtor(l,j,k,i),l=1,2)
+ enddo
+ enddo
+ enddo
+ endif
+ ELSE
+ do i=-nloctyp+1,nloctyp-1
+ do ii=1,3
+ do j=1,2
+ bnew1tor(ii,j,i)=bnew1(ii,j,i)
+ enddo
+ enddo
+ do ii=1,3
+ do j=1,2
+ bnew2tor(ii,j,i)=bnew2(ii,j,i)
+ enddo
+ enddo
+ do ii=1,3
+ ccnewtor(ii,1,i)=ccnew(ii,1,i)
+ ccnewtor(ii,2,i)=ccnew(ii,2,i)
+ ddnewtor(ii,1,i)=ddnew(ii,1,i)
+ ddnewtor(ii,2,i)=ddnew(ii,2,i)
+ enddo
+ enddo
+ ENDIF !SPLIT_FOURIER_TOR
+#else
+ allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
+ allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
+ allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
+ allocate(b(13,-nloctyp-1:nloctyp+1))
+ if (lprint) &
+ write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
+ do i=0,nloctyp-1
+ read (ifourier,*,end=115,err=115)
+ read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
+ if (lprint) then
+ write (iout,*) 'Type ',onelet(iloctyp(i))
+ write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
+ endif
+ if (i.gt.0) then
+ b(2,-i)= b(2,i)
+ b(3,-i)= b(3,i)
+ b(4,-i)=-b(4,i)
+ b(5,-i)=-b(5,i)
+ endif
+!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)
+!cc B1tilde(1,i) = b(3,i)
+!cc B1tilde(2,i) =-b(5,i)
+!c B1tilde(1,-i) =-b(3,i)
+!c B1tilde(2,-i) =b(5,i)
+!cc b1tilde(1,i)=0.0d0
+!cc b1tilde(2,i)=0.0d0
+!cc B2(1,i) = b(2,i)
+!cc B2(2,i) = b(4,i)
+!c B2(1,-i) =b(2,i)
+!c B2(2,-i) =-b(4,i)
+
+!c b2(1,i)=0.0d0
+!c b2(2,i)=0.0d0
+ CCold(1,1,i)= b(7,i)
+ CCold(2,2,i)=-b(7,i)
+ CCold(2,1,i)= b(9,i)
+ CCold(1,2,i)= b(9,i)
+ CCold(1,1,-i)= b(7,i)
+ CCold(2,2,-i)=-b(7,i)
+ CCold(2,1,-i)=-b(9,i)
+ CCold(1,2,-i)=-b(9,i)
+!c CC(1,1,i)=0.0d0
+!c CC(2,2,i)=0.0d0
+!c CC(2,1,i)=0.0d0
+!c CC(1,2,i)=0.0d0
+!c Ctilde(1,1,i)= CCold(1,1,i)
+!c Ctilde(1,2,i)= CCold(1,2,i)
+!c Ctilde(2,1,i)=-CCold(2,1,i)
+!c Ctilde(2,2,i)=-CCold(2,2,i)
+!c CC(1,1,i)=0.0d0
+!c CC(2,2,i)=0.0d0
+!c CC(2,1,i)=0.0d0
+!c CC(1,2,i)=0.0d0
+!c Ctilde(1,1,i)= CCold(1,1,i)
+!c Ctilde(1,2,i)= CCold(1,2,i)
+!c Ctilde(2,1,i)=-CCold(2,1,i)
+!c Ctilde(2,2,i)=-CCold(2,2,i)
+
+!c Ctilde(1,1,i)=0.0d0
+!c Ctilde(1,2,i)=0.0d0
+!c Ctilde(2,1,i)=0.0d0
+!c Ctilde(2,2,i)=0.0d0
+ DDold(1,1,i)= b(6,i)
+ DDold(2,2,i)=-b(6,i)
+ DDold(2,1,i)= b(8,i)
+ DDold(1,2,i)= b(8,i)
+ DDold(1,1,-i)= b(6,i)
+ DDold(2,2,-i)=-b(6,i)
+ DDold(2,1,-i)=-b(8,i)
+ DDold(1,2,-i)=-b(8,i)
+!c DD(1,1,i)=0.0d0
+!c DD(2,2,i)=0.0d0
+!c DD(2,1,i)=0.0d0
+!c DD(1,2,i)=0.0d0
+!c Dtilde(1,1,i)= DD(1,1,i)
+!c Dtilde(1,2,i)= DD(1,2,i)
+!c Dtilde(2,1,i)=-DD(2,1,i)
+!c Dtilde(2,2,i)=-DD(2,2,i)
+
+!c Dtilde(1,1,i)=0.0d0
+!c Dtilde(1,2,i)=0.0d0
+!c Dtilde(2,1,i)=0.0d0
+!c Dtilde(2,2,i)=0.0d0
+ EEold(1,1,i)= b(10,i)+b(11,i)
+ EEold(2,2,i)=-b(10,i)+b(11,i)
+ EEold(2,1,i)= b(12,i)-b(13,i)
+ EEold(1,2,i)= b(12,i)+b(13,i)
+ EEold(1,1,-i)= b(10,i)+b(11,i)
+ EEold(2,2,-i)=-b(10,i)+b(11,i)
+ EEold(2,1,-i)=-b(12,i)+b(13,i)
+ EEold(1,2,-i)=-b(12,i)-b(13,i)
+ write(iout,*) "TU DOCHODZE"
+ print *,"JESTEM"
+!c ee(1,1,i)=1.0d0
+!c ee(2,2,i)=1.0d0
+!c ee(2,1,i)=0.0d0
+!c ee(1,2,i)=0.0d0
+!c ee(2,1,i)=ee(1,2,i)
+ enddo
+ if (lprint) then
+ write (iout,*)
+ write (iout,*) &
+ "Coefficients of the cumulants (independent of valence angles)"
+ do i=-nloctyp+1,nloctyp-1
+ write (iout,*) 'Type ',onelet(iloctyp(i))
+ write (iout,*) 'B1'
+ write(iout,'(2f10.5)') B(3,i),B(5,i)
+ write (iout,*) 'B2'
+ write(iout,'(2f10.5)') B(2,i),B(4,i)
+ write (iout,*) 'CC'
+ do j=1,2
+ write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
+ enddo
+ write(iout,*) 'DD'
+ do j=1,2
+ write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
+ enddo
+ write(iout,*) 'EE'
+ do j=1,2
+ write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
+ enddo
+ enddo
+ endif
+#endif
+
+
#ifdef CRYST_TOR
!
! Read torsional parameters in old format
!
! Read torsional parameters
!
+ IF (TOR_MODE.eq.0) THEN
allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
read (itorp,*,end=113,err=113) ntortyp
!el from energy module---------
enddo
enddo
endif
+ ELSE IF (TOR_MODE.eq.1) THEN
+
+!C read valence-torsional parameters
+ read (itorp,*,end=121,err=121) ntortyp
+ nkcctyp=ntortyp
+ write (iout,*) "Valence-torsional parameters read in ntortyp",&
+ ntortyp
+ read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
+ write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
+#ifndef NEWCORR
+ do i=1,ntyp1
+ itype2loc(i)=itortyp(i)
+ enddo
+#endif
+ do i=-ntyp,-1
+ itortyp(i)=-itortyp(-i)
+ enddo
+ do i=-ntortyp+1,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+!C first we read the cos and sin gamma parameters
+ read (itorp,'(13x,a)',end=121,err=121) string
+ write (iout,*) i,j,string
+ read (itorp,*,end=121,err=121) &
+ nterm_kcc(j,i),nterm_kcc_Tb(j,i)
+!C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
+ do k=1,nterm_kcc(j,i)
+ do l=1,nterm_kcc_Tb(j,i)
+ do ll=1,nterm_kcc_Tb(j,i)
+ read (itorp,*,end=121,err=121) ii,jj,kk, &
+ v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ ELSE
+#ifdef NEWCORR
+!c AL 4/8/16: Calculate coefficients from one-body parameters
+ ntortyp=nloctyp
+ allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
+ allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
+ allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
+ allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
+ allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
+
+ do i=-ntyp1,ntyp1
+ print *,i,itortyp(i)
+ itortyp(i)=itype2loc(i)
+ enddo
+ write (iout,*) &
+ "Val-tor parameters calculated from cumulant coefficients ntortyp"&
+ ,ntortyp
+ do i=-ntortyp+1,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+ nterm_kcc(j,i)=2
+ nterm_kcc_Tb(j,i)=3
+ do k=1,nterm_kcc_Tb(j,i)
+ do l=1,nterm_kcc_Tb(j,i)
+ v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
+ +bnew1tor(k,2,i)*bnew2tor(l,2,j)
+ v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
+ +bnew1tor(k,2,i)*bnew2tor(l,1,j)
+ enddo
+ enddo
+ do k=1,nterm_kcc_Tb(j,i)
+ do l=1,nterm_kcc_Tb(j,i)
+#ifdef CORRCD
+ v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
+ -ccnewtor(k,2,i)*ddnewtor(l,2,j))
+ v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
+ +ccnewtor(k,1,i)*ddnewtor(l,2,j))
+#else
+ v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
+ -ccnewtor(k,2,i)*ddnewtor(l,2,j))
+ v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
+ +ccnewtor(k,1,i)*ddnewtor(l,2,j))
+#endif
+ enddo
+ enddo
+!c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
+ enddo
+ enddo
+#else
+ write (iout,*) "TOR_MODE>1 only with NEWCORR"
+ stop
+#endif
+ ENDIF ! TOR_MODE
+
+ if (tor_mode.gt.0 .and. lprint) then
+!c Print valence-torsional parameters
+ write (iout,'(a)') &
+ "Parameters of the valence-torsional potentials"
+ do i=-ntortyp+1,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+ write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
+ write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
+ do k=1,nterm_kcc(j,i)
+ do l=1,nterm_kcc_Tb(j,i)
+ do ll=1,nterm_kcc_Tb(j,i)
+ write (iout,'(3i5,2f15.4)')&
+ k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ endif
#endif
allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
! interaction energy of the Gly, Ala, and Pro prototypes.
!
-
- if (lprint) then
- write (iout,*)
- write (iout,*) "Coefficients of the cumulants"
- endif
- read (ifourier,*) nloctyp
-!write(iout,*) "nloctyp",nloctyp
-!el from module energy-------
- allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
- allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
- allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
- allocate(cc(2,2,-nloctyp-1:nloctyp+1))
- allocate(dd(2,2,-nloctyp-1:nloctyp+1))
- allocate(ee(2,2,-nloctyp-1:nloctyp+1))
- allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
- allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
-! el
- b1(1,:)=0.0d0
- b1(2,:)=0.0d0
-!--------------------------------
-
- do i=0,nloctyp-1
- read (ifourier,*,end=115,err=115)
- read (ifourier,*,end=115,err=115) (buse(ii),ii=1,13)
- if (lprint) then
- write (iout,*) 'Type',i
- write (iout,'(a,i2,a,f10.5)') ('buse(',ii,')=',buse(ii),ii=1,13)
- endif
- B1(1,i) = buse(3)
- B1(2,i) = buse(5)
- B1(1,-i) = buse(3)
- B1(2,-i) = -buse(5)
-! buse1(1,i)=0.0d0
-! buse1(2,i)=0.0d0
- B1tilde(1,i) = buse(3)
- B1tilde(2,i) =-buse(5)
- B1tilde(1,-i) =-buse(3)
- B1tilde(2,-i) =buse(5)
-! buse1tilde(1,i)=0.0d0
-! buse1tilde(2,i)=0.0d0
- B2(1,i) = buse(2)
- B2(2,i) = buse(4)
- B2(1,-i) =buse(2)
- B2(2,-i) =-buse(4)
-
-! buse2(1,i)=0.0d0
-! buse2(2,i)=0.0d0
- CC(1,1,i)= buse(7)
- CC(2,2,i)=-buse(7)
- CC(2,1,i)= buse(9)
- CC(1,2,i)= buse(9)
- CC(1,1,-i)= buse(7)
- CC(2,2,-i)=-buse(7)
- CC(2,1,-i)=-buse(9)
- CC(1,2,-i)=-buse(9)
-! CC(1,1,i)=0.0d0
-! CC(2,2,i)=0.0d0
-! CC(2,1,i)=0.0d0
-! CC(1,2,i)=0.0d0
- Ctilde(1,1,i)=buse(7)
- Ctilde(1,2,i)=buse(9)
- Ctilde(2,1,i)=-buse(9)
- Ctilde(2,2,i)=buse(7)
- Ctilde(1,1,-i)=buse(7)
- Ctilde(1,2,-i)=-buse(9)
- Ctilde(2,1,-i)=buse(9)
- Ctilde(2,2,-i)=buse(7)
-
-! Ctilde(1,1,i)=0.0d0
-! Ctilde(1,2,i)=0.0d0
-! Ctilde(2,1,i)=0.0d0
-! Ctilde(2,2,i)=0.0d0
- DD(1,1,i)= buse(6)
- DD(2,2,i)=-buse(6)
- DD(2,1,i)= buse(8)
- DD(1,2,i)= buse(8)
- DD(1,1,-i)= buse(6)
- DD(2,2,-i)=-buse(6)
- DD(2,1,-i)=-buse(8)
- DD(1,2,-i)=-buse(8)
-! DD(1,1,i)=0.0d0
-! DD(2,2,i)=0.0d0
-! DD(2,1,i)=0.0d0
-! DD(1,2,i)=0.0d0
- Dtilde(1,1,i)=buse(6)
- Dtilde(1,2,i)=buse(8)
- Dtilde(2,1,i)=-buse(8)
- Dtilde(2,2,i)=buse(6)
- Dtilde(1,1,-i)=buse(6)
- Dtilde(1,2,-i)=-buse(8)
- Dtilde(2,1,-i)=buse(8)
- Dtilde(2,2,-i)=buse(6)
-
-! Dtilde(1,1,i)=0.0d0
-! Dtilde(1,2,i)=0.0d0
-! Dtilde(2,1,i)=0.0d0
-! Dtilde(2,2,i)=0.0d0
- EE(1,1,i)= buse(10)+buse(11)
- EE(2,2,i)=-buse(10)+buse(11)
- EE(2,1,i)= buse(12)-buse(13)
- EE(1,2,i)= buse(12)+buse(13)
- EE(1,1,-i)= buse(10)+buse(11)
- EE(2,2,-i)=-buse(10)+buse(11)
- EE(2,1,-i)=-buse(12)+buse(13)
- EE(1,2,-i)=-buse(12)-buse(13)
-
-! ee(1,1,i)=1.0d0
-! ee(2,2,i)=1.0d0
-! ee(2,1,i)=0.0d0
-! ee(1,2,i)=0.0d0
-! ee(2,1,i)=ee(1,2,i)
- enddo
- if (lprint) then
- do i=1,nloctyp
- write (iout,*) 'Type',i
- write (iout,*) 'B1'
- write(iout,*) B1(1,i),B1(2,i)
- write (iout,*) 'B2'
- write(iout,*) B2(1,i),B2(2,i)
- write (iout,*) 'CC'
- do j=1,2
- write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
- enddo
- write(iout,*) 'DD'
- do j=1,2
- write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
- enddo
- write(iout,*) 'EE'
- do j=1,2
- write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
- enddo
- enddo
- endif
!
! Read electrostatic-interaction parameters
!
allocate(wstate(4,ntyp,ntyp))
allocate(dhead(2,2,ntyp,ntyp))
allocate(nstate(ntyp,ntyp))
+ allocate(debaykap(ntyp,ntyp))
+
if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+
do i=1,ntyp
do j=1,i
! write (*,*) "Im in ALAB", i, " ", j
epshead(i,j),sig0head(i,j), &
rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
alphapol(i,j),alphapol(j,i), &
- (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
+ (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j)
! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
END DO
END DO
wquad(i,j) = wquad(j,i)
epsintab(i,j) = epsintab(j,i)
+ debaykap(i,j)=debaykap(j,i)
! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
END DO
END DO
sigeps=dsign(1.0D0,epsij)
epsij=dabs(epsij)
aa_aq(i,j)=epsij*rrij*rrij
- print *,"ADASKO",epsij,rrij,expon
+! print *,"ADASKO",epsij,rrij,expon
bb_aq(i,j)=-sigeps*epsij*rrij
aa_aq(j,i)=aa_aq(i,j)
bb_aq(j,i)=bb_aq(i,j)
118 write (iout,*) "Error reading SCp interaction parameters."
goto 999
119 write (iout,*) "Error reading SCCOR parameters"
+ go to 999
+ 121 write (iout,*) "Error in Czybyshev parameters"
999 continue
#ifdef MPI
call MPI_Finalize(Ierror)
logical :: lprn=.true.,fail
real(kind=8),dimension(3) :: e1,e2,e3
real(kind=8) :: dcj,efree_temp
- character(len=3) :: seq,res
+ character(len=3) :: seq,res,res2
character(len=5) :: atom
character(len=80) :: card
real(kind=8),dimension(3,20) :: sccor
if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
! Fish out the ATOM cards.
! write(iout,*) 'card',card(1:20)
+! print *,"ATU ",card(1:6), CARD(3:6)
+! print *,card
if (index(card(1:4),'ATOM').gt.0) then
read (card(12:16),*) atom
! write (iout,*) "! ",atom," !",ires
enddo
else
call sccenter(ires_old,iii,sccor)
- endif
+ endif !unres_pdb
iii=0
- endif
+ endif !ind_pdb
! Start new residue.
if (res.eq.'Cl-' .or. res.eq.'Na+') then
ires=ires_old
ishift=ishift-1
itype(1,1)=ntyp1
nres_molec(molecule)=nres_molec(molecule)+1
- endif
+ endif ! Gly
ires=ires-ishift+ishift1
ires_old=ires
! write (iout,*) "ishift",ishift," ires",ires,&
ishift=ishift-(ires-ishift+ishift1-ires_old-1)
ires=ires-ishift+ishift1
ires_old=ires
- endif
+ endif ! Na Cl
! print *,'atom',ires,atom
if (res.eq.'ACE' .or. res.eq.'NHE') then
itype(ires,1)=10
! nres_molec(molecule)=nres_molec(molecule)+1
else
molecule=2
- itype(ires,molecule)=rescode(ires,res(2:3),0,molecule)
+ res2=res(2:3)
+ itype(ires,molecule)=rescode(ires,res2,0,molecule)
! nres_molec(molecule)=nres_molec(molecule)+1
read (card(19:19),'(a1)') sugar
isugar=sugarcode(sugar,ires)
! print *,"ires=",ires,istype(ires)
! endif
- endif
- endif
+ endif ! atom.eq.CA
+ endif !ACE
else
ires=ires-ishift+ishift1
- endif
+ endif !ires_old
! write (iout,*) "ires_old",ires_old," ires",ires
if (card(27:27).eq."A" .or. card(27:27).eq."B") then
! ishift1=ishift1+1
read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
endif
endif
+! print *,"IONS",ions,card(1:6)
else if ((ions).and.(card(1:6).eq.'HETATM')) then
if (firstion.eq.0) then
firstion=1
enddo
else
call sccenter(ires,iii,sccor)
- endif
- endif
+ endif ! unres_pdb
+ endif !firstion
read (card(12:16),*) atom
- print *,"HETATOM", atom
+! print *,"HETATOM", atom
read (card(18:20),'(a3)') res
if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
(atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
res=res(2:3)//' '
itype(ires,molecule)=rescode(ires,res,0,molecule)
read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
- endif
+ endif! NA
endif !atom
enddo
10 write (iout,'(a,i5)') ' Number of residues found: ',ires
! e2(3)=0.0d0
! endif
do j=1,3
- c(j,i)=c(j,i+1)-1.9d0*e2(j)
+! c(j,i)=c(j,i+1)-1.9d0*e2(j)
+ c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
enddo
! else !unres_pdb
do j=1,3
e2(3)=0.0d0
endif
do j=1,3
- c(j,1)=c(j,2)-1.9d0*e2(j)
+! c(j,1)=c(j,2)-1.9d0*e2(j)
+ c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
enddo
else
do j=1,3
call readi(controlcard,'SHIELD',shield_mode,0)
!C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
write(iout,*) "shield_mode",shield_mode
+ call readi(controlcard,'TORMODE',tor_mode,0)
+!C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write(iout,*) "torsional and valence angle mode",tor_mode
+
!C Varibles set size of box
with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
protein=index(controlcard,"PROTEIN").gt.0
projects / unres4.git / blobdiff