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)
dsc_inv(i)=1.0D0/dsc(i)
endif
enddo
+! ip(1)=0.0001d0
+! isc(:,1)=0.0001d0
#endif
read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
do i=1,ntyp_molec(2)
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
+! ip(5)=0.2
! isc(5)=0.2
read (iion,*) ncatprotparm
allocate(catprm(ncatprotparm,4))
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
enddo
endif
enddo
+#ifdef SC_END
+ allocate(nterm_scend(2,ntyp))
+ allocate(arotam_end(0:6,2,ntyp))
+ nterm_scend=0
+ arotam_end=0.0d0
+ read (irotam_end,*) ijunk
+!c write (iout,*) "ijunk",ijunk
+ do i=1,ntyp
+ if (i.eq.10) cycle
+ do j=1,2
+ read (irotam_end,'(a)')
+ read (irotam_end,*) nterm_scend(j,i)
+!c write (iout,*) "i",i," j",j," nterm",nterm_scend(j,i)
+ do k=0,nterm_scend(j,i)
+ read (irotam_end,*) ijunk,arotam_end(k,j,i)
+!c write (iout,*) "k",k," arotam",arotam_end(k,j,i)
+ enddo
+ enddo
+ enddo
+!c lprint=.true.
+ if (lprint) then
+ write (iout,'(a)') &
+ "Parameters of the local potentials of sidechain ends"
+ do i=1,ntyp
+ write (iout,'(5x,9x,2hp-,a3,6x,9x,a3,2h-p)')&
+ restyp(i,1),restyp(i,1)
+ do j=0,max0(nterm_scend(1,i),nterm_scend(2,i))
+ write (iout,'(i5,2f20.10)') &
+ j,arotam_end(j,1,i),arotam_end(j,2,i)
+ enddo
+ enddo
+ endif
+!c lprint=.false.
+#endif
+
!---------reading nucleic acid parameters for rotamers-------------------
allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
do i=1,ntyp_molec(2)
#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
sigiso2(i,j)=sigiso2(j,i)
! print *,"ATU",sigma(j,i),sigma(i,j),i,j
nstate(i,j) = nstate(j,i)
- dtail(1,i,j) = dtail(1,j,i)
- dtail(2,i,j) = dtail(2,j,i)
+ dtail(1,i,j) = dtail(2,j,i)
+ dtail(2,i,j) = dtail(1,j,i)
DO k = 1, 4
alphasur(k,i,j) = alphasur(k,j,i)
wstate(k,i,j) = wstate(k,j,i)
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
-
- if(me.eq.king) then
- write (iout,'(/a)') "Disulfide bridge parameters:"
- write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
- write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
- write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
- write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
- ' v3ss:',v3ss
+! 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)=72.0d0
+ mp(4)=0.0d0
+ isc(i,4)=0.d0
+ enddo
+ ip(4)=0.0
+ msc(ntyp_molec(4)+1,4)=0.1d0
+ !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
- if (shield_mode.gt.0) then
- pi=4.0D0*datan(1.0D0)
-!C VSolvSphere the volume of solving sphere
- print *,pi,"pi"
+ 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,-1:ntyp_molec(5)),epsheadcat(ntyp,-1:ntyp_molec(5)),sig0headcat(ntyp,-1:ntyp_molec(5)))
+ allocate(alphapolcat2(ntyp,-1:ntyp_molec(5)))
+ allocate(sigiso1cat(ntyp,-1:ntyp_molec(5)),rborn1cat(ntyp,-1:ntyp_molec(5)),rborn2cat(ntyp,-1:ntyp_molec(5)),sigmap1cat(ntyp,-1:ntyp_molec(5)))
+ allocate(sigmap2cat(ntyp,-1:ntyp_molec(5)),sigiso2cat(ntyp,-1:ntyp_molec(5)))
+ allocate(chis1cat(ntyp,-1:ntyp_molec(5)),chis2cat(ntyp,-1:ntyp_molec(5)),wquadcat(ntyp,-1:ntyp_molec(5)),chipp1cat(ntyp,-1:ntyp_molec(5)),chipp2cat(ntyp,-1:ntyp_molec(5)))
+ allocate(epsintabcat(ntyp,-1:ntyp_molec(5)))
+ allocate(dtailcat(2,ntyp,-1:ntyp_molec(5)))
+ allocate(alphasurcat(4,ntyp,-1:ntyp_molec(5)),alphisocat(4,ntyp,-1:ntyp_molec(5)))
+ allocate(wqdipcat(2,ntyp,-1:ntyp_molec(5)))
+ allocate(wstatecat(4,ntyp,-1:ntyp_molec(5)))
+ allocate(dheadcat(2,2,ntyp,-1:ntyp_molec(5)))
+ allocate(nstatecat(ntyp,-1:ntyp_molec(5)))
+ allocate(debaykapcat(ntyp,-1:ntyp_molec(5)))
+
+ if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,-1:ntyp1))
+ if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,-1:ntyp1))
+! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
+ if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,-1:ntyp1)) !(ntyp,ntyp)
+ if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,-1: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
+ if (j.eq.0) cycle
+ 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), & !6
+
+ (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&!12
+! chiscat(i,j),chiscat(j,i), &
+ chis1cat(i,j),chis2cat(i,j), &
+
+ nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !19 !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),&!23
+ dtailcat(1,i,j),dtailcat(2,i,j), &
+ epsheadcat(i,j),sig0headcat(i,j), &!27
+!wdipcat = w1 , w2
+! rborncat(i,j),rborncat(j,i),&
+ rborn1cat(i,j),rborn2cat(i,j),&
+ (wqdipcat(k,i,j),k=1,2), &!31
+ alphapolcat(i,j),alphapolcat2(j,i), &!33
+ (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,-1:ntyp_molec(5)),&
+ bb_aq_cat(-ntyp:ntyp,-1:ntyp_molec(5)))
+ do i=1,ntyp
+ do j=-1,ntyp_molec(5)
+ if (j.eq.0) cycle
+ 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
+ write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
+ write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
+ write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
+ ' v3ss:',v3ss
+ endif
+
+!------------MARTINI-PROTEIN-parameters-------------------------
+ allocate(alphapolmart(ntyp,ntyp),epsheadmart(ntyp,ntyp_molec(4)),sig0headmart(ntyp,ntyp_molec(4)))
+ allocate(alphapolmart2(ntyp,ntyp))
+ allocate(sigiso1mart(ntyp,ntyp_molec(4)),rborn1mart(ntyp,ntyp_molec(4)),rborn2mart(ntyp,ntyp_molec(4)),sigmap1mart(ntyp,ntyp_molec(4)))
+ allocate(sigmap2mart(ntyp,ntyp_molec(4)),sigiso2mart(ntyp,ntyp_molec(4)))
+ allocate(chis1mart(ntyp,ntyp_molec(4)),chis2mart(ntyp,ntyp_molec(4)),wquadmart(ntyp,ntyp_molec(4)),chipp1mart(ntyp,ntyp_molec(4)),chipp2mart(ntyp,ntyp_molec(4)))
+ allocate(epsintabmart(ntyp,ntyp_molec(4)))
+ allocate(dtailmart(2,ntyp,ntyp_molec(4)))
+ allocate(alphasurmart(4,ntyp,ntyp_molec(4)),alphisomart(4,ntyp,ntyp_molec(4)))
+ allocate(wqdipmart(2,ntyp,ntyp_molec(4)))
+ allocate(wstatemart(4,ntyp,ntyp_molec(4)))
+ allocate(dheadmart(2,2,ntyp,ntyp_molec(4)))
+ allocate(nstatemart(ntyp,ntyp_molec(4)))
+ allocate(debaykapmart(ntyp,ntyp_molec(4)))
+
+ if (.not.allocated(epsmart)) allocate (epsmart(0:ntyp1,ntyp1))
+ if (.not.allocated(sigmamart)) allocate(sigmamart(0:ntyp1,ntyp1))
+! if (.not.allocated(chimart)) allomarte(chimart(ntyp1,ntyp1)) !(ntyp,ntyp)
+ if (.not.allocated(chi1mart)) allocate(chi1mart(ntyp1,ntyp1)) !(ntyp,ntyp)
+ if (.not.allocated(chi2mart)) allocate(chi2mart(ntyp1,ntyp1)) !(ntyp,ntyp)
+
+!DIR$ NOUNROLL
+ do i=1,ntyp-3 ! there are phosporylated missing
+ do j=1,ntyp_molec(4) ! this is without Zn will be modified for ALL tranistion metals
+! do j=1,ntyp_molec(5)
+ print *,"lipmart",i,j
+! write (*,*) "Im in ALAB", i, " ", j
+ read(imartprot,*) &
+ epsmart(i,j),sigmamart(i,j), &
+! chimart(i,j),chimart(j,i),chippmart(i,j),chippmart(j,i), &
+ chi1mart(i,j),chi2mart(i,j),chipp1mart(i,j),chipp2mart(i,j), & !6
+
+ (alphasurmart(k,i,j),k=1,4),sigmap1mart(i,j),sigmap2mart(i,j),&!12
+! chismart(i,j),chismart(j,i), &
+ chis1mart(i,j),chis2mart(i,j), &
+
+ nstatemart(i,j),(wstatemart(k,i,j),k=1,4), & !19 !5 w tej lini - 1 integer pierwszy
+ dheadmart(1,1,i,j),dheadmart(1,2,i,j),dheadmart(2,1,i,j),dheadmart(2,2,i,j),&!23
+ dtailmart(1,i,j),dtailmart(2,i,j), &
+ epsheadmart(i,j),sig0headmart(i,j), &!27
+!wdipmart = w1 , w2
+! rbornmart(i,j),rbornmart(j,i),&
+ rborn1mart(i,j),rborn2mart(i,j),&
+ (wqdipmart(k,i,j),k=1,2), &!31
+ alphapolmart(i,j),alphapolmart2(j,i), &!33
+ (alphisomart(k,i,j),k=1,4),sigiso1mart(i,j),sigiso2mart(i,j),epsintabmart(i,j),debaykapmart(i,j)
+ enddo
+ enddo
+ allocate(aa_aq_mart(-ntyp:ntyp,ntyp_molec(4)),&
+ bb_aq_mart(-ntyp:ntyp,ntyp_molec(4)))
+ do i=1,ntyp-3 ! still no phophorylated residues
+ do j=1,ntyp_molec(4)
+ if (j.eq.0) cycle
+ epsij=epsmart(i,j)
+ rrij=sigmamart(i,j)
+ rrij=rrij**expon
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa_aq_mart(i,j)=epsij*rrij*rrij
+ bb_aq_mart(i,j)=-sigeps*epsij*rrij
+ enddo
+ enddo
+
+
+
+
+
+
+
+
+
+
+
+ if (shield_mode.gt.0) then
+ pi=4.0D0*datan(1.0D0)
+!C VSolvSphere the volume of solving sphere
+ print *,pi,"pi"
!C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
!C there will be no distinction between proline peptide group and normal peptide
!C group in case of shielding parameters
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)
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
+ 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
! 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
+ 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)
cou=1
write (iout,*) "symetr", symetr
do i=1,nres
- lll=lll+1
+ lll=lll+1
! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
- if (i.gt.1) then
- if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
- chain_length=lll-1
- kkk=kkk+1
+! if (i.gt.1) then
+! if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
+! chain_length=lll-1
+! kkk=kkk+1
! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
- lll=1
- endif
- endif
+! lll=1
+! endif
+! endif
do j=1,3
cref(j,i,cou)=c(j,i)
cref(j,i+nres,cou)=c(j,i+nres)
endif
enddo
enddo
- write (iout,*) chain_length
- if (chain_length.eq.0) chain_length=nres
- do j=1,3
- chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
- chain_rep(j,chain_length+nres,symetr) &
- =chain_rep(j,chain_length+nres,1)
- enddo
! diagnostic
! write (iout,*) "spraw lancuchy",chain_length,symetr
! do i=1,4
dc(j,0)=c(j,1)
enddo
- if (symetr.gt.1) then
- call permut(symetr)
- nperm=1
- do i=1,symetr
- nperm=nperm*i
- enddo
- do i=1,nperm
- write(iout,*) (tabperm(i,kkk),kkk=1,4)
- enddo
- do i=1,nperm
- cou=0
- do kkk=1,symetr
- icha=tabperm(i,kkk)
- write (iout,*) i,icha
- do lll=1,chain_length
- cou=cou+1
- if (cou.le.nres) then
- do j=1,3
- kupa=mod(lll,chain_length)
- iprzes=(kkk-1)*chain_length+lll
- if (kupa.eq.0) kupa=chain_length
- write (iout,*) "kupa", kupa
- cref(j,iprzes,i)=chain_rep(j,kupa,icha)
- cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
- enddo
- endif
- enddo
- enddo
- enddo
- endif
+! if (symetr.gt.1) then
+! call permut(symetr)
+! nperm=1
+! do i=1,symetr
+! nperm=nperm*i
+! enddo
+! do i=1,nperm
+! write(iout,*) (tabperm(i,kkk),kkk=1,4)
+! enddo
+! do i=1,nperm
+! cou=0
+! do kkk=1,symetr
+! icha=tabperm(i,kkk)
+! write (iout,*) i,icha
+! do lll=1,chain_length
+! cou=cou+1
+! if (cou.le.nres) then
+! do j=1,3
+! kupa=mod(lll,chain_length)
+! iprzes=(kkk-1)*chain_length+lll
+! if (kupa.eq.0) kupa=chain_length
+! write (iout,*) "kupa", kupa
+! cref(j,iprzes,i)=chain_rep(j,kupa,icha)
+! cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
+! enddo
+! endif
+! enddo
+! enddo
+! enddo
+! endif
!-koniec robienia kopii
! diag
do kkk=1,nperm
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
+
+ write(iout,*) "R_CUT_ELE=",r_cut_ele,rlamb_ele
+ 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)
+
+ call reada(controlcard,"DELTA",graddelta,1.0d-5)
+! 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
! print *,"Processor",myrank," opened file ITHEP"
call getenv_loc('ROTPAR',rotname)
open (irotam,file=rotname,status='old',action='read')
+#ifdef SC_END
+ call getenv_loc('ROTPAR_END',rotname_end)
+ open (irotam_end,file=rotname_end,status='old',action='read')
+#endif
! print *,"Processor",myrank," opened file IROTAM"
call getenv_loc('TORPAR',torname)
open (itorp,file=torname,status='old',action='read')
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 (ithep,file=thetname,status='old')
call getenv_loc('ROTPAR',rotname)
open (irotam,file=rotname,status='old')
+#ifdef SC_END
+ call getenv_loc('ROTPAR_END',rotname_end)
+ open (irotam_end,file=rotname_end,status='old')
+#endif
call getenv_loc('TORPAR',torname)
open (itorp,file=torname,status='old')
call getenv_loc('TORDPAR',tordname)
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)
open (ithep,file=thetname,status='old',action='read')
call getenv_loc('ROTPAR',rotname)
open (irotam,file=rotname,status='old',action='read')
+#ifdef SC_END
+ call getenv_loc('ROTPAR_END',rotname_end)
+ open (irotam_end,file=rotname_end,status='old',action='read')
+#endif
call getenv_loc('TORPAR',torname)
open (itorp,file=torname,status='old',action='read')
call getenv_loc('TORDPAR',tordname)
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('LIPPROT',lipprotname)
+ open (imartprot,file=lipprotname,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
projects / unres4.git / blobdiff