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
+ sigmasctube,krad2,ract
+ integer :: ichir1,ichir2,ijunk,irdiff
+ character*3 string
+ character*80 temp1,mychar
! 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))
msc(:,:)=0.0d0
isc(:,:)=0.0d0
- allocate(msc(ntyp+1,5)) !(ntyp+1)
- allocate(isc(ntyp+1,5)) !(ntyp+1)
- allocate(restok(ntyp+1,5)) !(ntyp+1)
+ allocate(msc(-ntyp-1:ntyp+1,5)) !(ntyp+1)
+ allocate(isc(-ntyp-1:ntyp+1,5)) !(ntyp+1)
+ allocate(restok(-ntyp-1:ntyp+1,5)) !(ntyp+1)
read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
do i=1,ntyp_molec(1)
enddo
enddo
endif
+
+
+
+ if (.not.allocated(ichargecat)) &
+ allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
+ ichargecat(:)=0
+ if (oldion.eq.1) then
do i=1,ntyp_molec(5)
- read(iion,*) msc(i,5),restok(i,5)
+ read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
print *,msc(i,5),restok(i,5)
enddo
ip(5)=0.2
read (iion,*) (catprm(i,k),i=1,ncatprotparm)
enddo
print *, catprm
+ endif
+ allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
+ do i=1,ntyp_molec(5)
+ do j=1,ntyp_molec(2)
+ write(iout,*) i,j
+ read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
+ enddo
+ enddo
+ write(*,'(3(5x,a6)11(7x,a6))') "w1 ","w2 ","epslj ","pis1 ", &
+ "sigma0","epsi0 ","chi1 ","chip1 ","sig ","b1 ","b2 ", &
+ "b3 ","b4 ","chis1 "
+ do i=1,ntyp_molec(5)
+ do j=1,ntyp_molec(2)
+ write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
+ restyp(i,5),"-",restyp(j,2)
+ enddo
+ enddo
! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
!----------------------------------------------------
allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
print *,liptranene(i)
enddo
close(iliptranpar)
+! water parmaters entalphy
+ allocate(awaterenta(0:400))
+ allocate(bwaterenta(0:400))
+ allocate(cwaterenta(0:400))
+ allocate(dwaterenta(0:400))
+ allocate(awaterentro(0:400))
+ allocate(bwaterentro(0:400))
+ allocate(cwaterentro(0:400))
+ allocate(dwaterentro(0:400))
+
+ read(iwaterwater,*) mychar
+ read(iwaterwater,*) ract,awaterenta(0),bwaterenta(0),&
+ cwaterenta(0),dwaterenta(0)
+ do i=1,398
+ read(iwaterwater,*) ract,awaterenta(i),bwaterenta(i),&
+ cwaterenta(i),dwaterenta(i)
+ irdiff=int((ract-2.06d0)*50.0d0)+1
+ if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff
+ enddo
+! water parmaters entrophy
+ read(iwaterwater,*) mychar
+ read(iwaterwater,*) ract,awaterentro(0),bwaterentro(0),&
+ cwaterentro(0),dwaterentro(0)
+ do i=1,398
+ read(iwaterwater,*) ract,awaterentro(i),bwaterentro(i),&
+ cwaterentro(i),dwaterentro(i)
+ irdiff=int((ract-2.06d0)*50.0d0)+1
+ if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff
+ enddo
+
#ifdef CRYST_THETA
!
! 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))
+ bnew1=0.0d0
+ bnew2=0.0d0
+ ccnew=0.0d0
+ ddnew=0.0d0
+ e0new=0.0d0
+ eenew=0.0d0
+ bnew1tor=0.0d0
+ bnew2tor=0.0d0
+ ccnewtor=0.0d0
+ ddnewtor=0.0d0
+ e0newtor=0.0d0
+ eenewtor=0.0d0
+ 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
+#ifndef NEWCORR
+ do i=1,ntyp1
+ itype2loc(i)=itortyp(i)
+ 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
si=-1.0d0
do k=1,nterm_sccor(i,j)
+ print *,"test",i,j,k,l
read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
v2sccor(k,l,i,j)
v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
! 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
read(isidep,*) &
- eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
- (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
- chis(i,j),chis(j,i), &
- nstate(i,j),(wstate(k,i,j),k=1,4), &
- dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
- dtail(1,i,j),dtail(2,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)
-! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
+ eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
+ (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
+ chis(i,j),chis(j,i), & !2 w tej linii
+ nstate(i,j),(wstate(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
+ dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& ! 4 w tej linii
+ dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
+ epshead(i,j),sig0head(i,j), & ! 2 w tej linii
+ rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & ! 5 w tej linii
+ alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
+ (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
+ IF ((LaTeX).and.(i.gt.24)) then
+ write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
+ eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
+ (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
+ chis(i,j),chis(j,i) !2 w tej linii
+ endif
+ print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i), wqdip(1,i,j)
+ END DO
+ END DO
+ do i=1,ntyp
+ read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
+ enddo
+ do i=1,ntyp
+ do j=1,i
+ IF ((LaTeX).and.(i.gt.24)) then
+ write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
+ dhead(1,1,i,j),dhead(2,1,i,j),& ! 2 w tej linii
+ dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
+ epshead(i,j),sig0head(i,j), & ! 2 w tej linii
+ rborn(i,j),rborn(j,i), & ! 3 w tej linii
+ alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
+ (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
+ endif
END DO
END DO
DO i = 1, ntyp
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
endif
+
+
! goto 50
!--------------------- GBV potential -----------------------------------
case(5)
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)
allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
-
+ write (iout,*) "ESCBASEPARM"
do i=1,ntyp_molec(1)
- do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
- write (*,*) "Im in ", i, " ", j
+ do j=1,ntyp_molec(2) ! without U then we will take T for U
+! write (*,*) "Im in ", i, " ", j
read(isidep_scbase,*) &
eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
- write(*,*) "eps",eps_scbase(i,j)
+! write(*,*) "eps",eps_scbase(i,j)
read(isidep_scbase,*) &
(alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
chis_scbase(i,j,1),chis_scbase(i,j,2)
read(isidep_scbase,*) &
alphapol_scbase(i,j), &
epsintab_scbase(i,j)
+ if (chi_scbase(i,j,2).gt.0.9) chi_scbase(i,j,2)=0.9
+ if (chi_scbase(i,j,1).gt.0.9) chi_scbase(i,j,1)=0.9
+ if (chipp_scbase(i,j,2).gt.0.9) chipp_scbase(i,j,2)=0.9
+ if (chipp_scbase(i,j,1).gt.0.9) chipp_scbase(i,j,1)=0.9
+ if (chi_scbase(i,j,2).lt.-0.9) chi_scbase(i,j,2)=-0.9
+ if (chi_scbase(i,j,1).lt.-0.9) chi_scbase(i,j,1)=-0.9
+ if (chipp_scbase(i,j,2).lt.-0.9) chipp_scbase(i,j,2)=-0.9
+ if (chipp_scbase(i,j,1).lt.-0.9) chipp_scbase(i,j,1)=-0.9
+ write(iout,*) &
+ eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
+ chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
+ write(*,*) "eps",eps_scbase(i,j)
+ write(iout,*) &
+ (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
+ chis_scbase(i,j,1),chis_scbase(i,j,2)
+ write(iout,*) &
+ dhead_scbasei(i,j), &
+ dhead_scbasej(i,j), &
+ rborn_scbasei(i,j),rborn_scbasej(i,j)
+ write(iout,*) &
+ (wdipdip_scbase(k,i,j),k=1,3), &
+ (wqdip_scbase(k,i,j),k=1,2)
+ write(iout,*) &
+ alphapol_scbase(i,j), &
+ epsintab_scbase(i,j)
+
END DO
+ j=4
+ write(iout,*) &
+ eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
+ chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
+ write(*,*) "eps",eps_scbase(i,j)
+ write(iout,*) &
+ (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
+ chis_scbase(i,j,1),chis_scbase(i,j,2)
+ write(iout,*) &
+ dhead_scbasei(i,j), &
+ dhead_scbasej(i,j), &
+ rborn_scbasei(i,j),rborn_scbasej(i,j)
+ write(iout,*) &
+ (wdipdip_scbase(k,i,j),k=1,3), &
+ (wqdip_scbase(k,i,j),k=1,2)
+ write(iout,*) &
+ alphapol_scbase(i,j), &
+ epsintab_scbase(i,j)
+
END DO
allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
allocate(chis_pepbase(ntyp_molec(2),2))
allocate(wdipdip_pepbase(3,ntyp_molec(2)))
+ write (iout,*) "EPEPBASEPARM"
- do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
+ do j=1,ntyp_molec(2) ! without U then we will take T for U
write (*,*) "Im in ", i, " ", j
read(isidep_pepbase,*) &
eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
+ if (chi_pepbase(j,2).gt.0.9) chi_pepbase(j,2)=0.9
+ if (chi_pepbase(j,1).gt.0.9) chi_pepbase(j,1)=0.9
+ if (chipp_pepbase(j,2).gt.0.9) chipp_pepbase(j,2)=0.9
+ if (chipp_pepbase(j,1).gt.0.9) chipp_pepbase(j,1)=0.9
+ if (chi_pepbase(j,2).lt.-0.9) chi_pepbase(j,2)=-0.9
+ if (chi_pepbase(j,1).lt.-0.9) chi_pepbase(j,1)=-0.9
+ if (chipp_pepbase(j,2).lt.-0.9) chipp_pepbase(j,2)=-0.9
+ if (chipp_pepbase(j,1).lt.-0.9) chipp_pepbase(j,1)=-0.9
+
write(*,*) "eps",eps_pepbase(j)
read(isidep_pepbase,*) &
(alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
chis_pepbase(j,1),chis_pepbase(j,2)
read(isidep_pepbase,*) &
(wdipdip_pepbase(k,j),k=1,3)
+ write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
+ chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
+ write(iout,*) &
+ (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
+ chis_pepbase(j,1),chis_pepbase(j,2)
+ write(iout,*) &
+ (wdipdip_pepbase(k,j),k=1,3)
+
END DO
+ j=4
+ write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
+ chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
+ write(iout,*) &
+ (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
+ chis_pepbase(j,1),chis_pepbase(j,2)
+ write(iout,*) &
+ (wdipdip_pepbase(k,j),k=1,3)
+
allocate(aa_pepbase(ntyp_molec(2)))
allocate(bb_pepbase(ntyp_molec(2)))
read(isidep_scpho,*) &
epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
alphi_scpho(j)
+ if (chi_scpho(j,2).gt.0.9) chi_scpho(j,2)=0.9
+ if (chi_scpho(j,1).gt.0.9) chi_scpho(j,1)=0.9
+ if (chipp_scpho(j,2).gt.0.9) chipp_scpho(j,2)=0.9
+ if (chipp_scpho(j,1).gt.0.9) chipp_scpho(j,1)=0.9
+ if (chi_scpho(j,2).lt.-0.9) chi_scpho(j,2)=-0.9
+ if (chi_scpho(j,1).lt.-0.9) chi_scpho(j,1)=-0.9
+ if (chipp_scpho(j,2).lt.-0.9) chipp_scpho(j,2)=-0.9
+ if (chipp_scpho(j,1).lt.-0.9) chipp_scpho(j,1)=-0.9
+
END DO
allocate(aa_scpho(ntyp_molec(1)))
! v1ss=0.0d0
! v2ss=0.0d0
! v3ss=0.0d0
-
+! MARTINI PARAMETER
+ allocate(ichargelipid(ntyp_molec(4)))
+ allocate(lip_angle_force(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
+ allocate(lip_angle_angle(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
+ allocate(lip_bond(ntyp_molec(4),ntyp_molec(4)))
+ allocate(lip_eps(ntyp_molec(4),ntyp_molec(4)))
+ allocate(lip_sig(ntyp_molec(4),ntyp_molec(4)))
+ kjtokcal=0.2390057361
+ krad=57.295779513
+ !HERE THE MASS of MARTINI
+ write(*,*) "before MARTINI PARAM"
+ do i=1,ntyp_molec(4)
+ msc(i,4)=0.0d0
+ mp(4)=72.0d0
+ isc(i,4)=0.d0
+ enddo
+ ip(4)=0.0
+ !relative dielectric constant = 15 for implicit screening
+ k_coulomb_lip=332.0d0/15.0d0
+ !kbond = 1250 kJ/(mol*nm*2)
+ kbondlip=1250.0d0*kjtokcal/100.0d0
+ krad2=krad**2.0
+ lip_angle_force=0.0d0
+ if (DRY_MARTINI.gt.0) then
+ lip_angle_force(3,12,12)=35.0*kjtokcal!*krad2
+ lip_angle_force(3,12,18)=35.0*kjtokcal!*krad2
+ lip_angle_force(3,18,16)=35.0*kjtokcal!*krad2
+ lip_angle_force(12,18,16)=35.0*kjtokcal!*krad2
+ lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
+ lip_angle_force(16,18,18)=35.0*kjtokcal!*krad2
+ lip_angle_force(12,18,18)=35.0*kjtokcal!*krad2
+ lip_angle_force(18,18,18)=35.0*kjtokcal!*krad2
+ else
+ lip_angle_force(3,12,12)=25.0*kjtokcal!*krad2
+ lip_angle_force(3,12,18)=25.0*kjtokcal!*krad2
+ lip_angle_force(3,18,16)=25.0*kjtokcal!*krad2
+ lip_angle_force(12,18,16)=25.0*kjtokcal!*krad2
+ lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
+ lip_angle_force(16,18,18)=25.0*kjtokcal!*krad2
+ lip_angle_force(12,18,18)=25.0*kjtokcal!*krad2
+ lip_angle_force(18,18,18)=25.0*kjtokcal!*krad2
+ endif
+ lip_angle_angle=0.0d0
+ lip_angle_angle(3,12,12)=120.0/krad
+ lip_angle_angle(3,12,18)=180.0/krad
+ lip_angle_angle(3,18,16)=180.0/krad
+ lip_angle_angle(12,18,16)=180.0/krad
+ lip_angle_angle(18,16,18)=120.0/krad
+ lip_angle_angle(16,18,18)=180.0/krad
+ lip_angle_angle(12,18,18)=180.0/krad
+ lip_angle_angle(18,18,18)=180.0/krad
+ read(ilipbond,*) temp1
+ do i=1,18
+ read(ilipbond,*) temp1, lip_bond(i,1), &
+ lip_bond(i,2),lip_bond(i,3),lip_bond(i,4),lip_bond(i,5), &
+ lip_bond(i,6),lip_bond(i,7),lip_bond(i,8),lip_bond(i,9), &
+ lip_bond(i,10),lip_bond(i,11),lip_bond(i,12),lip_bond(i,13), &
+ lip_bond(i,14),lip_bond(i,15),lip_bond(i,16),lip_bond(i,17), &
+ lip_bond(i,18)
+ do j=1,18
+ lip_bond(i,j)=lip_bond(i,j)*10
+ enddo
+ enddo
+
+ read(ilipnonbond,*) (ichargelipid(i),i=1,ntyp_molec(4))
+ read(ilipnonbond,*) temp1
+ do i=1,18
+ read(ilipnonbond,*) temp1, lip_eps(i,1), &
+ lip_eps(i,2),lip_eps(i,3),lip_eps(i,4),lip_eps(i,5), &
+ lip_eps(i,6),lip_eps(i,7),lip_eps(i,8),lip_eps(i,9), &
+ lip_eps(i,10),lip_eps(i,11),lip_eps(i,12),lip_eps(i,13), &
+ lip_eps(i,14),lip_eps(i,15),lip_eps(i,16),lip_eps(i,17), &
+ lip_eps(i,18)
+! write(*,*) i, lip_eps(i,18)
+ do j=1,18
+ lip_eps(i,j)=lip_eps(i,j)*kjtokcal
+ enddo
+ enddo
+ read(ilipnonbond,*) temp1
+ do i=1,18
+ read(ilipnonbond,*) temp1,lip_sig(i,1), &
+ lip_sig(i,2),lip_sig(i,3),lip_sig(i,4),lip_sig(i,5), &
+ lip_sig(i,6),lip_sig(i,7),lip_sig(i,8),lip_sig(i,9), &
+ lip_sig(i,10),lip_sig(i,11),lip_sig(i,12),lip_sig(i,13), &
+ lip_sig(i,14),lip_sig(i,15),lip_sig(i,16),lip_sig(i,17), &
+ lip_sig(i,18)
+ do j=1,18
+ lip_sig(i,j)=lip_sig(i,j)*10.0
+ enddo
+ enddo
+ write(*,*) "after MARTINI PARAM"
+
+! Ions by Aga
+
+ allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
+ allocate(alphapolcat2(ntyp,ntyp))
+ allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
+ allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
+ allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
+ allocate(epsintabcat(ntyp,ntyp))
+ allocate(dtailcat(2,ntyp,ntyp))
+ allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
+ allocate(wqdipcat(2,ntyp,ntyp))
+ allocate(wstatecat(4,ntyp,ntyp))
+ allocate(dheadcat(2,2,ntyp,ntyp))
+ allocate(nstatecat(ntyp,ntyp))
+ allocate(debaykapcat(ntyp,ntyp))
+
+ if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
+ if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
+! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
+ if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
+ if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
+
+
+ if (.not.allocated(ichargecat))&
+ allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
+ write(*,*) "before ions",oldion
+ ichargecat(:)=0
+
+! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
+ if (oldion.eq.0) then
+ if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
+ allocate(icharge(1:ntyp1))
+ read(iion,*) (icharge(i),i=1,ntyp)
+ else
+ read(iion,*) ijunk
+ endif
+ print *,ntyp_molec(5)
+ do i=-ntyp_molec(5),ntyp_molec(5)
+ read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
+ print *,msc(i,5),restok(i,5)
+ enddo
+ ip(5)=0.2
+ ! mp(5)=0.2
+ pstok(5)=3.0
+!DIR$ NOUNROLL
+ do j=1,ntyp_molec(5)-1 ! this is without Zn will be modified for ALL tranistion metals
+ do i=1,ntyp
+! do j=1,ntyp_molec(5)
+! write (*,*) "Im in ALAB", i, " ", j
+ read(iion,*) &
+ epscat(i,j),sigmacat(i,j), &
+! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
+ chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
+
+ (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
+! chiscat(i,j),chiscat(j,i), &
+ chis1cat(i,j),chis2cat(i,j), &
+
+ nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
+ dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
+ dtailcat(1,i,j),dtailcat(2,i,j), &
+ epsheadcat(i,j),sig0headcat(i,j), &
+!wdipcat = w1 , w2
+! rborncat(i,j),rborncat(j,i),&
+ rborn1cat(i,j),rborn2cat(i,j),&
+ (wqdipcat(k,i,j),k=1,2), &
+ alphapolcat(i,j),alphapolcat2(j,i), &
+ (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
+
+ if (chi1cat(i,j).gt.0.9) write (*,*) "WTF ANISO", i,j, chi1cat(i,j)
+! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
+! if (i.eq.1) then
+! write (iout,*) 'i= ', i, ' j= ', j
+! write (iout,*) 'epsi0= ', epscat(1,j)
+! write (iout,*) 'sigma0= ', sigmacat(1,j)
+! write (iout,*) 'chi(i,j)= ', chicat(1,j)
+! write (iout,*) 'chi(j,i)= ', chicat(j,1)
+! write (iout,*) 'chip(i,j)= ', chippcat(1,j)
+! write (iout,*) 'chip(j,i)= ', chippcat(j,1)
+! write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
+! write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
+! write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
+! write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
+! write (iout,*) 'sig1= ', sigmap1cat(1,j)
+! write (iout,*) 'chis(i,j)= ', chiscat(1,j)
+! write (iout,*) 'chis(j,i)= ', chiscat(j,1)
+! write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
+! write (iout,*) 'a1= ', rborncat(j,1)
+! write (iout,*) 'a2= ', rborncat(1,j)
+! write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
+! write (iout,*) 'alphapol1= ', alphapolcat(1,j)
+! write (iout,*) 'w1= ', wqdipcat(1,1,j)
+! write (iout,*) 'w2= ', wqdipcat(2,1,j)
+! endif
+
+!
+! If ((i.eq.1).and.(j.eq.27)) then
+! write (iout,*) 'SEP'
+! Write (iout,*) 'w1= ', wqdipcat(1,1,27)
+! Write (iout,*) 'w2= ', wqdipcat(2,1,27)
+! endif
+
+ END DO
+ END DO
+ allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
+ do i=1,ntyp
+ do j=1,ntyp_molec(5)
+ epsij=epscat(i,j)
+ rrij=sigmacat(i,j)
+ rrij=rrij**expon
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa_aq_cat(i,j)=epsij*rrij*rrij
+ bb_aq_cat(i,j)=-sigeps*epsij*rrij
+ enddo
+ enddo
+
+ do i=1,ntyp
+ do j=1,ntyp_molec(5)-1
+ if (i.eq.10) then
+ write (iout,*) 'i= ', i, ' j= ', j
+ write (iout,*) 'epsi0= ', epscat(i,j)
+ write (iout,*) 'sigma0= ', sigmacat(i,j)
+ write (iout,*) 'chi1= ', chi1cat(i,j)
+ write (iout,*) 'chi1= ', chi2cat(i,j)
+ write (iout,*) 'chip1= ', chipp1cat(i,j)
+ write (iout,*) 'chip2= ', chipp2cat(i,j)
+ write (iout,*) 'alphasur1= ', alphasurcat(1,i,j)
+ write (iout,*) 'alphasur2= ', alphasurcat(2,i,j)
+ write (iout,*) 'alphasur3= ', alphasurcat(3,i,j)
+ write (iout,*) 'alphasur4= ', alphasurcat(4,i,j)
+ write (iout,*) 'sig1= ', sigmap1cat(i,j)
+ write (iout,*) 'sig2= ', sigmap2cat(i,j)
+ write (iout,*) 'chis1= ', chis1cat(i,j)
+ write (iout,*) 'chis1= ', chis2cat(i,j)
+ write (iout,*) 'nstatecat(i,j)= ', nstatecat(i,j)
+ write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,i,j)
+ write (iout,*) 'dhead= ', dheadcat(1,1,i,j)
+ write (iout,*) 'dhead2= ', dheadcat(1,2,i,j)
+ write (iout,*) 'a1= ', rborn1cat(i,j)
+ write (iout,*) 'a2= ', rborn2cat(i,j)
+ write (iout,*) 'epsin= ', epsintabcat(i,j), epsintabcat(j,i)
+ write (iout,*) 'alphapol1= ', alphapolcat(i,j)
+ write (iout,*) 'alphapol2= ', alphapolcat2(i,j)
+ write (iout,*) 'w1= ', wqdipcat(1,i,j)
+ write (iout,*) 'w2= ', wqdipcat(2,i,j)
+ write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(i,j)
+ endif
+
+ If ((i.eq.1).and.(j.eq.27)) then
+ write (iout,*) 'SEP'
+ Write (iout,*) 'w1= ', wqdipcat(1,1,27)
+ Write (iout,*) 'w2= ', wqdipcat(2,1,27)
+ endif
+
+ enddo
+ enddo
+
+ endif
+! read number of Zn2+
+! here two denotes the Zn2+ and Cu2+
+ write(iout,*) "before TRANPARM"
+ allocate(aomicattr(0:3,2))
+ allocate(athetacattran(0:6,5,2))
+ allocate(agamacattran(3,5,2))
+ allocate(acatshiftdsc(5,2))
+ allocate(bcatshiftdsc(5,2))
+ allocate(demorsecat(5,2))
+ allocate(alphamorsecat(5,2))
+ allocate(x0catleft(5,2))
+ allocate(x0catright(5,2))
+ allocate(x0cattrans(5,2))
+ allocate(ntrantyp(2))
+ do i=1,1 ! currently only Zn2+
+
+ read(iiontran,*) ntrantyp(i)
+!ntrantyp=4
+!| ao0 ao1 ao2 ao3
+!ASP| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
+!CYS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0left x0right x0transi
+!GLU| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -0.5
+!HIS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
+ read(iiontran,*) (aomicattr(j,i),j=0,3)
+ do j=1,ntrantyp(i)
+ read (iiontran,*,err=123,end=123) (agamacattran(k,j,i),k=1,3),&
+ (athetacattran(k,j,i),k=0,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),&
+ demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),&
+ x0cattrans(j,i)
+ enddo
+ enddo
if(me.eq.king) then
write (iout,'(/a)') "Disulfide bridge parameters:"
write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
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"
+ go to 999
+ 123 write(iout,*) "Error in transition metal parameters"
999 continue
#ifdef MPI
call MPI_Finalize(Ierror)
character(len=3) :: seq,res,res2
character(len=5) :: atom
character(len=80) :: card
- real(kind=8),dimension(3,20) :: sccor
+ real(kind=8),dimension(3,40) :: sccor
integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
integer :: isugar,molecprev,firstion
character*1 :: sugar
itype(ires_old,molecule)=ntyp1_molec(molecule)
itype(ires_old-1,molecule)=ntyp1_molec(molecule)
nres_molec(molecule)=nres_molec(molecule)+2
+! if (molecule.eq.4) ires=ires+2
ibeg=2
! write (iout,*) "Chain ended",ires,ishift,ires_old
if (unres_pdb) then
! iii=0
endif
iii=0
+ else if (card(:3).eq.'BRA') then
+ molecule=4
+ ires=ires+1
+ ires_old=ires+1
+ itype(ires,molecule)=ntyp1_molec(molecule)-1
+ nres_molec(molecule)=nres_molec(molecule)+1
+
endif
! Read free energy
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
! if (atom.eq.'CA' .or. atom.eq.'CH3') then
+ if (card(14:16).eq.'LIP') then
+! reading lipid
+ if (ibeg.eq.1) then
+ molecule=4
+ ires=ires+1
+ nres_molec(molecule)=nres_molec(molecule)+1
+ itype(ires,molecule)=ntyp1_molec(molecule)
+ ibeg=0
+ endif
+ if (ibeg.eq.2) then
+ ibeg=0
+ ires=ires+2
+ endif
+
+ molecule=4
+ nres_molec(molecule)=nres_molec(molecule)+1
+ read (card(18:20),'(a3)') res
+ ires=ires+1
+ read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+ if (UNRES_PDB) then!
+ itype(ires,molecule)=rescode(ires,res,0,molecule)
+ else
+ itype(ires,molecule)=rescode_lip(res,ires)
+ endif
+ else
read (card(23:26),*) ires
read (card(18:20),'(a3)') res
! write (iout,*) "ires",ires,ires-ishift+ishift1,
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
! 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
+ endif !LIP
+! 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(1:2)
read (card(18:20),'(a3)') res
+ if (atom(3:3).eq.'H') cycle
if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
(atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
- .or.(atom(1:2).eq.'K ')) &
- then
+ .or.(atom(1:2).eq.'K ').or.(atom(1:2).eq.'ZN').or.&
+ (atom(1:2).eq.'O ')) then
ires=ires+1
+ print *,"I have water"
if (molecule.ne.5) molecprev=molecule
molecule=5
nres_molec(molecule)=nres_molec(molecule)+1
print *,"HERE",nres_molec(molecule)
- res=res(2:3)//' '
+ if (res.ne.'WAT') 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
! print *,"molecule",molecule
if ((itype(nres,1).ne.10)) then
nres=nres+1
+ nsup=nsup+1
if (molecule.eq.5) molecule=molecprev
itype(nres,molecule)=ntyp1_molec(molecule)
nres_molec(molecule)=nres_molec(molecule)+1
! enddo
! else
do j=1,3
- dcj=(c(j,nres-2)-c(j,nres-3))/2.0
- c(j,nres)=c(j,nres-1)+dcj
+ dcj=(c(j,nres-nres_molec(5)-2)-c(j,nres-nres_molec(5)-3))/2.0
+ c(j,nres)=c(j,nres-nres_molec(5)-1)+dcj
c(j,2*nres)=c(j,nres)
enddo
! endif
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
do j=1,3
c_temporary(j,seqalingbegin)=c(j,i)
c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
-
enddo
itype_temporary(seqalingbegin,k)=itype(i,k)
print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
do i=1,2*nres
do j=1,3
c(j,i)=c_temporary(j,i)
+ if (i.gt.nres) then
+ if ((molnum(i-nres).eq.4)) c(j,i)=0.0d0
+ endif
enddo
enddo
do k=1,5
istype(i)=istype_temp(i)
enddo
enddo
+ if ((itype(1,1).eq.ntyp1).and.itype(2,5).ne.0) then
+! I have only ions now dummy atoms in the system
+ molnum(1)=5
+ itype(1,5)=itype(2,5)
+ itype(1,1)=0
+ do j=1,3
+ c(j,1)=c(j,2)
+ enddo
+ do i=2,nres-1
+ itype(i,5)=itype(i+1,5)
+ do j=1,3
+ c(j,i)=c(j,i+1)
+ enddo
+ enddo
+ itype(nres,5)=0
+ nres=nres-1
+ nres_molec(1)=nres_molec(1)-1
+ endif
! if (itype(1,1).eq.ntyp1) then
! nsup=nsup-1
! nstart_sup=2
enddo
endif
-! print *,seqalingbegin,nres
+ print *,seqalingbegin,nres
if(.not.allocated(vbld)) then
allocate(vbld(2*nres))
do i=1,2*nres
allocate(dc_norm(3,0:2*nres+2))
dc_norm(:,:)=0.d0
endif
-
+ write(iout,*) "before int_from_cart"
call int_from_cart(.true.,.false.)
call sc_loc_geom(.false.)
+ write(iout,*) "after int_from_cart"
+
+
do i=1,nres
thetaref(i)=theta(i)
phiref(i)=phi(i)
enddo
+ write(iout,*) "after thetaref"
! do i=1,2*nres
! vbld_inv(i)=0.d0
! vbld(i)=0.d0
! enddo
! enddo
!
+! do i=1,2*nres
+! do j=1,3
+! chomo(j,i,k)=c(j,i)
+! enddo
+! enddo
+! write(iout,*) "after chomo"
+
if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
integer i
-
+ usampl=.false. ! the default value of usample should be 0
+! write(iout,*) "KURWA2",usampl
nglob_csa=0
eglob_csa=1d99
nmin_csa=0
call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
call reada(controlcard,'DRMS',drms,0.1D0)
+ call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
+ read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
+ out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
+ out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
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
ions=index(controlcard,"IONS").gt.0
+ fodson=index(controlcard,"FODSON").gt.0
+ call readi(controlcard,'OLDION',oldion,1)
nucleic=index(controlcard,"NUCLEIC").gt.0
write (iout,*) "with_theta_constr ",with_theta_constr
AFMlog=(index(controlcard,'AFM'))
! elemode = 0 is orignal UNRES electrostatics
! elemode = 1 is "Momo" potentials in progress
! elemode = 2 is in development EVALD
+
+
write (iout,*) TUBEmode,"TUBEMODE"
if (TUBEmode.gt.0) then
call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
call reada(controlcard,"RTUBE",tubeR0,0.0d0)
+ call reada(controlcard,"VNANO",velnanoconst,0.0d0)
call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
endif
! CUTOFFF ON ELECTROSTATICS
- call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
+ call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
write(iout,*) "R_CUT_ELE=",r_cut_ele
-! Lipidic parameters
+ call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
+ call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
+ call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
+
+! Lipidec parameters
call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
if (lipthick.gt.0.0d0) then
! enddo
buff_shield=1.0d0
endif
+ itime_mat=0
return
end subroutine read_control
!-----------------------------------------------------------------------------
!
! Read MD settings
!
- use control_data, only: r_cut,rlamb,out1file
+ use control_data, only: r_cut,rlamb,out1file,r_cut_mart,rlamb_mart
use energy_data
use geometry_data, only: pi
use MPI_data
call readi(controlcard,"NSTEP",n_timestep,1000000)
call readi(controlcard,"NTWE",ntwe,100)
call readi(controlcard,"NTWX",ntwx,1000)
+ call readi(controlcard,"NFOD",nfodstep,100)
call reada(controlcard,"DT",d_time,1.0d-1)
call reada(controlcard,"DVMAX",dvmax,2.0d1)
call reada(controlcard,"DAMAX",damax,1.0d1)
rest = index(controlcard,"REST").gt.0
tbf = index(controlcard,"TBF").gt.0
usampl = index(controlcard,"USAMPL").gt.0
+! write(iout,*) "KURWA",usampl
mdpdb = index(controlcard,"MDPDB").gt.0
call reada(controlcard,"T_BATH",t_bath,300.0d0)
call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
print_compon = index(controlcard,"PRINT_COMPON").gt.0
rattle = index(controlcard,"RATTLE").gt.0
preminim=(index(controlcard,'PREMINIM').gt.0)
+ forceminim=(index(controlcard,'FORCEMINIM').gt.0)
write (iout,*) "PREMINIM ",preminim
dccart=(index(controlcard,'CART').gt.0)
if (preminim) call read_minim
enddo
if(me.eq.king.or..not.out1file)then
+ do j=1,5
write (iout,'(/2a/)') &
"Radii of site types and friction coefficients and std's of",&
" stochastic forces of fully exposed sites"
- write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
+ write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(j),stdfp*dsqrt(gamp(j))
+
do i=1,ntyp
- write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
- gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
+ write (iout,'(a5,f5.2,2f10.5)') restyp(i,j),restok(i,j),&
+ gamsc(i,j),stdfsc(i,j)*dsqrt(gamsc(i,j))
enddo
+ enddo
endif
else if (tbf) then
if(me.eq.king.or..not.out1file)then
open (itube,file=tubename,status='old',action='read')
call getenv_loc('IONPAR',ionname)
open (iion,file=ionname,status='old',action='read')
-
+ call getenv_loc('IONPAR_TRAN',iontranname)
+ open (iiontran,file=iontranname,status='old',action='read')
! print *,"Processor",myrank," opened file ISIDEP"
! print *,"Processor",myrank," opened parameter files"
#elif (defined G77)
open (itube,file=tubename,status='old')
call getenv_loc('IONPAR',ionname)
open (iion,file=ionname,status='old')
+ call getenv_loc('IONPAR_NUCL',ionnuclname)
+ open (iionnucl,file=ionnuclname,status='old')
+ call getenv_loc('IONPAR_TRAN',iontranname)
+ open (iiontran,file=iontranname,status='old')
+ call getenv_loc('WATWAT',iwaterwatername)
+ open (iwaterwater,file=iwaterwatername,status='old')
+ call getenv_loc('WATPROT',iwaterscname)
+ open (iwatersc,file=iwaterscname,status='old')
+
#else
open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
readonly)
call getenv_loc('LIPTRANPAR',liptranname)
open (iliptranpar,file=liptranname,status='old',action='read')
+ call getenv_loc('LIPBOND',lipbondname)
+ open (ilipbond,file=lipbondname,status='old',action='read')
+ call getenv_loc('LIPNONBOND',lipnonbondname)
+ open (ilipnonbond,file=lipnonbondname,status='old',action='read')
call getenv_loc('TUBEPAR',tubename)
open (itube,file=tubename,status='old',action='read')
call getenv_loc('IONPAR',ionname)
open (iion,file=ionname,status='old',action='read')
+ call getenv_loc('IONPAR_NUCL',ionnuclname)
+ open (iionnucl,file=ionnuclname,status='old',action='read')
+ call getenv_loc('IONPAR_TRAN',iontranname)
+ open (iiontran,file=iontranname,status='old',action='read')
+ call getenv_loc('WATWAT',iwaterwatername)
+ open (iwaterwater,file=iwaterwatername,status='old',action='read')
+ call getenv_loc('WATPROT',iwaterscname)
+ open (iwatersc,file=iwaterscname,status='old',action='read')
#ifndef CRYST_SC
call getenv_loc('ROTPARPDB',rotname_pdb)
end subroutine write_stat_thread
!-----------------------------------------------------------------------------
#endif
+ subroutine readpdb_template(k)
+! Read the PDB file for read_constr_homology with read2sigma
+! and convert the peptide geometry into virtual-chain geometry.
+! implicit none
+! include 'DIMENSIONS'
+! include 'COMMON.LOCAL'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.NAMES'
+! include 'COMMON.CONTROL'
+! include 'COMMON.FRAG'
+! include 'COMMON.SETUP'
+ use compare_data, only:nhfrag,nbfrag
+ integer :: i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity, &
+ ishift_pdb,ires_ca
+ logical lprn /.false./,fail
+ real(kind=8), dimension (3):: e1,e2,e3
+ real(kind=8) :: dcj,efree_temp
+ character*3 seq,res
+ character*5 atom
+ character*80 card
+ real(kind=8), dimension (3,40) :: sccor
+! integer rescode
+ integer, dimension (:), allocatable :: iterter
+ if(.not.allocated(iterter))allocate(iterter(nres))
+ do i=1,nres
+ iterter(i)=0
+ enddo
+ ibeg=1
+ ishift1=0
+ ishift=0
+ write (2,*) "UNRES_PDB",unres_pdb
+ ires=0
+ ires_old=0
+ iii=0
+ lsecondary=.false.
+ nhfrag=0
+ nbfrag=0
+ do
+ read (ipdbin,'(a80)',end=10) card
+ if (card(:3).eq.'END') then
+ goto 10
+ else if (card(:3).eq.'TER') then
+! End current chain
+ ires_old=ires+2
+ itype(ires_old-1,1)=ntyp1
+ iterter(ires_old-1)=1
+#if defined(WHAM_RUN) || defined(CLUSTER)
+ if (ires_old.lt.nres) then
+#endif
+ itype(ires_old,1)=ntyp1
+ iterter(ires_old)=1
+#if defined(WHAM_RUN) || defined(CLUSTER)
+ endif
+#endif
+ ibeg=2
+! write (iout,*) "Chain ended",ires,ishift,ires_old
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires,iii,sccor)
+ endif
+ endif
+! Fish out the ATOM cards.
+ if (index(card(1:4),'ATOM').gt.0) then
+ read (card(12:16),*) atom
+! write (iout,*) "! ",atom," !",ires
+! if (atom.eq.'CA' .or. atom.eq.'CH3') then
+ read (card(23:26),*) ires
+ read (card(18:20),'(a3)') res
+! write (iout,*) "ires",ires,ires-ishift+ishift1,
+! & " ires_old",ires_old
+! write (iout,*) "ishift",ishift," ishift1",ishift1
+! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+ if (ires-ishift+ishift1.ne.ires_old) then
+! Calculate the CM of the preceding residue.
+ if (ibeg.eq.0) then
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires_old)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires_old,iii,sccor)
+ endif
+ iii=0
+ endif
+! Start new residue.
+ if (res.eq.'Cl-' .or. res.eq.'Na+') then
+ ires=ires_old
+ cycle
+ else if (ibeg.eq.1) then
+! write (iout,*) "BEG ires",ires
+ ishift=ires-1
+ if (res.ne.'GLY' .and. res.ne. 'ACE') then
+ ishift=ishift-1
+ itype(1,1)=ntyp1
+ endif
+ ires=ires-ishift+ishift1
+ ires_old=ires
+! write (iout,*) "ishift",ishift," ires",ires,
+! & " ires_old",ires_old
+! write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
+ ibeg=0
+ else if (ibeg.eq.2) then
+! Start a new chain
+ ishift=-ires_old+ires-1
+ ires=ires_old+1
+! write (iout,*) "New chain started",ires,ishift
+ ibeg=0
+ else
+ ishift=ishift-(ires-ishift+ishift1-ires_old-1)
+ ires=ires-ishift+ishift1
+ ires_old=ires
+ endif
+ if (res.eq.'ACE' .or. res.eq.'NHE') then
+ itype(ires,1)=10
+ else
+ itype(ires,1)=rescode(ires,res,0,1)
+ endif
+ else
+ ires=ires-ishift+ishift1
+ endif
+! write (iout,*) "ires_old",ires_old," ires",ires
+! if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+! ishift1=ishift1+1
+! endif
+! write (2,*) "ires",ires," res ",res," ity",ity
+ if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
+ res.eq.'NHE'.and.atom(:2).eq.'HN') then
+ read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+! write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
+#ifdef DEBUG
+ write (iout,'(2i3,2x,a,3f8.3)') &
+ ires,itype(ires,1),res,(c(j,ires),j=1,3)
+#endif
+ iii=iii+1
+ do j=1,3
+ sccor(j,iii)=c(j,ires)
+ enddo
+ if (ishift.ne.0) then
+ ires_ca=ires+ishift-ishift1
+ else
+ ires_ca=ires
+ endif
+! write (*,*) card(23:27),ires,itype(ires)
+ else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.&
+ atom.ne.'N' .and. atom.ne.'C' .and.&
+ atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.&
+ atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+! write (iout,*) "sidechain ",atom
+ iii=iii+1
+ read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+ endif
+ endif
+ enddo
+ 10 if(me.eq.king.or..not.out1file) &
+ write (iout,'(a,i5)') ' Nres: ',ires
+! Calculate dummy residue coordinates inside the "chain" of a multichain
+! system
+ nres=ires
+ write(2,*) "tutaj",ires,nres
+ do i=2,nres-1
+! write (iout,*) i,itype(i),itype(i+1)
+ if (itype(i,1).eq.ntyp1.and.iterter(i).eq.1) then
+ if (itype(i+1,1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
+! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
+! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
+! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
+ if (unres_pdb) then
+! 2/15/2013 by Adam: corrected insertion of the last dummy residue
+ call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif !fail
+ do j=1,3
+ c(j,i)=c(j,i-1)-1.9d0*e2(j)
+ enddo
+ else !unres_pdb
+ do j=1,3
+ dcj=(c(j,i-2)-c(j,i-3))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,i)=c(j,i-1)+dcj
+ c(j,nres+i)=c(j,i)
+ enddo
+ endif !unres_pdb
+ else !itype(i+1).eq.ntyp1
+ if (unres_pdb) then
+! 2/15/2013 by Adam: corrected insertion of the first dummy residue
+ call refsys(i+1,i+2,i+3,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,i)=c(j,i+1)-1.9d0*e2(j)
+ enddo
+ else !unres_pdb
+ do j=1,3
+ dcj=(c(j,i+3)-c(j,i+2))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,i)=c(j,i+1)-dcj
+ c(j,nres+i)=c(j,i)
+ enddo
+ endif !unres_pdb
+ endif !itype(i+1).eq.ntyp1
+ endif !itype.eq.ntyp1
+ enddo
+! Calculate the CM of the last side chain.
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires,iii,sccor)
+ endif
+ nsup=nres
+ nstart_sup=1
+ if (itype(nres,1).ne.10) then
+ nres=nres+1
+ itype(nres,1)=ntyp1
+ if (unres_pdb) then
+! 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)-1.9d0*e2(j)
+ enddo
+ else
+ do j=1,3
+ dcj=(c(j,nres-2)-c(j,nres-3))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,nres)=c(j,nres-1)+dcj
+ c(j,2*nres)=c(j,nres)
+ enddo
+ endif
+ endif
+ do i=2,nres-1
+ do j=1,3
+ c(j,i+nres)=dc(j,i)
+ enddo
+ enddo
+ do j=1,3
+ c(j,nres+1)=c(j,1)
+ c(j,2*nres)=c(j,nres)
+ 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
+! Copy the coordinates to reference coordinates
+! do i=1,2*nres
+! do j=1,3
+! cref(j,i)=c(j,i)
+! enddo
+! enddo
+! Calculate internal coordinates.
+ if (out_template_coord) then
+ write (iout,'(/a)') &
+ "Cartesian coordinates of the reference structure"
+ write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
+ "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+ do ires=1,nres
+ write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)')&
+ restyp(itype(ires,1),1),ires,(c(j,ires),j=1,3),&
+ (c(j,ires+nres),j=1,3)
+ enddo
+ endif
+! Calculate internal coordinates.
+ call int_from_cart(.true.,out_template_coord)
+ call sc_loc_geom(.false.)
+ do i=1,nres
+ thetaref(i)=theta(i)
+ phiref(i)=phi(i)
+ enddo
+ do i=1,nres-1
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ enddo
+ enddo
+ do i=2,nres-1
+ do j=1,3
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ enddo
+! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
+! & vbld_inv(i+nres)
+ enddo
+ do i=1,nres
+ do j=1,3
+ cref(j,i,1)=c(j,i)
+ cref(j,i+nres,1)=c(j,i+nres)
+ enddo
+ enddo
+ do i=1,2*nres
+ do j=1,3
+ chomo(j,i,k)=c(j,i)
+ enddo
+ enddo
+
+ return
+ end subroutine readpdb_template
+!-----------------------------------------------------------------------------
+!#endif
!-----------------------------------------------------------------------------
end module io_config