From: Bartlomiej Zaborowski Date: Wed, 23 Jan 2013 11:39:13 +0000 (+0100) Subject: Mergowanie adasko do bartek X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=7fcc0efe632e299b8b944facdca6b9320c8c78da;hp=-c;p=unres.git Mergowanie adasko do bartek Merge branch 'adasko' into bartek Conflicts: PARAM/pot_tor_G631_DIL_ext.parm bin/wham/wham_multparm-ham_rep-oldparm source/unres/src_MD-M/cinfo.f source/unres/src_MD/COMMON.SCCOR source/unres/src_MD/Makefile source/unres/src_MD/cinfo.f source/unres/src_MD/parmread.F source/wham/src-M/cinfo.f source/wham/src/Makefile source/wham/src/cinfo.f source/wham/src/enecalc1.F source/wham/src/energy_p_new.F source/wham/src/include_unres/COMMON.SCCOR source/wham/src/int_from_cart.f source/wham/src/parmread.F --- 7fcc0efe632e299b8b944facdca6b9320c8c78da diff --combined bin/wham/wham_multparm-ham_rep-oldparm index 6f60701,3e7d413..0000000 deleted file mode 100755,100755 Binary files differ diff --combined source/cluster/wham/src/readrtns.F index 8e63ff8,4ccad2c..c61569b --- a/source/cluster/wham/src/readrtns.F +++ b/source/cluster/wham/src/readrtns.F @@@ -125,7 -125,6 +125,7 @@@ C Read weights of the subsequent energ call reada(weightcard,'SCALSCP',scalscp,1.0d0) call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) + call reada(weightcard,'WSCCOR',wsccor,1.0D0) if (index(weightcard,'SOFT').gt.0) ipot=6 C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) @@@ -150,7 -149,7 +150,7 @@@ weights(18)=scal14 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3, - & wturn4,wturn6 + & wturn4,wturn6,wsccor 10 format (/'Energy-term weights (unscaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ @@@ -168,8 -167,7 +168,8 @@@ & 'WCORR6= ',f10.6,' (multi-body 6th order)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ - & 'WTURN6= ',f10.6,' (turns, 6th order)') + & 'WTURN6= ',f10.6,' (turns, 6th order)'/ + & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)') if (wcorr4.gt.0.0d0) then write (iout,'(/2a/)') 'Local-electrostatic type correlation ', & 'between contact pairs of peptide groups' @@@ -210,15 -208,15 +210,15 @@@ C Convert sequence to numeric cod do i=1,nres #ifdef PROCOR - if (itype(i).eq.21 .or. itype(i+1).eq.21) then + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then #else - if (itype(i).eq.21) then + if (itype(i).eq.ntyp1) then #endif itel(i)=0 #ifdef PROCOR - else if (itype(i+1).ne.20) then + else if (iabs(itype(i+1)).ne.20) then #else - else if (itype(i).ne.20) then + else if (iabs(itype(i)).ne.20) then #endif itel(i)=1 else @@@ -254,8 -252,8 +254,8 @@@ nnt=1 nct=nres print *,'NNT=',NNT,' NCT=',NCT - if (itype(1).eq.21) nnt=2 - if (itype(nres).eq.21) nct=nct-1 + if (itype(1).eq.ntyp1) nnt=2 + if (itype(nres).eq.ntyp1) nct=nct-1 if (nstart.lt.nnt) nstart=nnt if (nend.gt.nct .or. nend.eq.0) nend=nct write (iout,*) "nstart",nstart," nend",nend diff --combined source/unres/src_MD/energy_p_new_barrier.F index 9edadf8,c0ceab6..bb2c70b --- a/source/unres/src_MD/energy_p_new_barrier.F +++ b/source/unres/src_MD/energy_p_new_barrier.F @@@ -1137,8 -1137,8 +1137,8 @@@ c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) - itypi1=iabs(itype(i+1)) + itypi=itype(i) + itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@@ -1151,7 -1151,7 +1151,7 @@@ cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) do j=istart(i,iint),iend(i,iint) - itypj=iabs(itype(j)) + itypj=itype(j) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@@ -1314,8 -1314,8 +1314,8 @@@ c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) - itypi1=iabs(itype(i+1)) + itypi=itype(i) + itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@@ -1324,7 -1324,7 +1324,7 @@@ C Calculate SC interaction energy C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) - itypj=iabs(itype(j)) + itypj=itype(j) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@@ -1431,8 -1431,8 +1431,8 @@@ c els c endif ind=0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) - itypi1=iabs(itype(i+1)) + itypi=itype(i) + itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@@ -1447,7 -1447,7 +1447,7 @@@ do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) chi1=chi(itypi,itypj) @@@ -1567,8 -1567,8 +1567,8 @@@ c print *,'Entering EGB nnt=',nnt, c if (icall.eq.0) lprn=.false. ind=0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) - itypi1=iabs(itype(i+1)) + itypi=itype(i) + itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@@ -1585,7 -1585,7 +1585,7 @@@ do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=iabs(itype(j)) + itypj=itype(j) c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, @@@ -1726,8 -1726,8 +1726,8 @@@ c print *,'Entering EGB nnt=',nnt, c if (icall.eq.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) - itypi1=iabs(itype(i+1)) + itypi=itype(i) + itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@@ -1742,7 -1742,7 +1742,7 @@@ do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=iabs(itype(j)) + itypj=itype(j) c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) @@@ -2046,8 -2046,8 +2046,8 @@@ cd print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) - itypi1=iabs(itype(i+1)) + itypi=itype(i) + itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@@ -2058,7 -2058,7 +2058,7 @@@ cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) do j=istart(i,iint),iend(i,iint) - itypj=iabs(itype(j)) + itypj=itype(j) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@@ -4059,7 -4059,7 +4059,7 @@@ cd write (iout,*) 'iatscp_s=',iatscp do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) - itypj=iabs(itype(j)) + itypj=itype(j) C Uncomment following three lines for SC-p interactions c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi @@@ -4153,7 -4153,7 +4153,7 @@@ cd write (iout,*) 'iatscp_s=',iatscp do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) - itypj=iabs(itype(j)) + itypj=itype(j) C Uncomment following three lines for SC-p interactions c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi @@@ -4270,7 -4270,8 +4270,7 @@@ c write (iout,*) "i",i," ii",ii, c & dhpb(i),dhpb1(i),forcon(i) C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. - if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. iabs(itype(jjj - &)).eq.1) then + if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij cd write (iout,*) "eij",eij @@@ -4374,7 -4375,7 +4374,7 @@@ include 'COMMON.VAR' include 'COMMON.IOUNITS' double precision erij(3),dcosom1(3),dcosom2(3),gg(3) - itypi=iabs(itype(i)) + itypi=itype(i) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@@ -4383,7 -4384,7 +4383,7 @@@ dzi=dc_norm(3,nres+i) c dsci_inv=dsc_inv(itypi) dsci_inv=vbld_inv(nres+i) - itypj=iabs(itype(j)) + itypj=itype(j) c dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(nres+j) xj=c(1,nres+j)-xi @@@ -4477,7 -4478,7 +4477,7 @@@ c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included c do i=ibond_start,ibond_end - iti=iabs(itype(i)) + iti=itype(i) if (iti.ne.10) then nbi=nbondterm(iti) if (nbi.eq.1) then @@@ -4552,7 -4553,19 +4552,7 @@@ c write (*,'(a,i2)') 'EBEND ICG=',i do i=ithet_start,ithet_end C Zero the energy function and its derivative at 0 or pi. call splinthet(theta(i),0.5d0*delta,ss,ssd) - it=(itype(i-1)) - ichir1=isign(1,itype(i-2)) - ichir2=isign(1,itype(i)) - if (itype(i-2).eq.10) ichir1=isign(1,itype(i-1)) - if (itype(i).eq.10) ichir2=isign(1,itype(i-1)) - if (itype(i-1).eq.10) then - itype1=isign(10,itype(i-2)) - ichir11=isign(1,itype(i-2)) - ichir12=isign(1,itype(i-2)) - itype2=isign(10,itype(i)) - ichir21=isign(1,itype(i)) - ichir22=isign(1,itype(i)) - endif + it=itype(i-1) if (i.gt.3) then #ifdef OSF phii=phi(i) @@@ -4586,15 -4599,27 +4586,15 @@@ C dependent on the adjacent virtual-bon C In following comments this theta will be referred to as t_c. thet_pred_mean=0.0d0 do k=1,2 - athetk=athet(k,it,ichir1,ichir2) - bthetk=bthet(k,it,ichir1,ichir2) - if (it.eq.10) then - athetk=athet(k,itype1,ichir11,ichir12) - bthetk=bthet(k,itype2,ichir21,ichir22) - endif + athetk=athet(k,it) + bthetk=bthet(k,it) thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k) enddo dthett=thet_pred_mean*ssd thet_pred_mean=thet_pred_mean*ss+a0thet(it) C Derivatives of the "mean" values in gamma1 and gamma2. - dthetg1=(-athet(1,it,ichir1,ichir2)*y(2) - &+athet(2,it,ichir1,ichir2)*y(1))*ss - dthetg2=(-bthet(1,it,ichir1,ichir2)*z(2) - & +bthet(2,it,ichir1,ichir2)*z(1))*ss - if (it.eq.10) then - dthetg1=(-athet(1,itype1,ichir11,ichir12)*y(2) - &+athet(2,itype1,ichir11,ichir12)*y(1))*ss - dthetg2=(-bthet(1,itype2,ichir21,ichir22)*z(2) - & +bthet(2,itype2,ichir21,ichir22)*z(1))*ss - endif + dthetg1=(-athet(1,it)*y(2)+athet(2,it)*y(1))*ss + dthetg2=(-bthet(1,it)*z(2)+bthet(2,it)*z(1))*ss if (theta(i).gt.pi-delta) then call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0, & E_tc0) @@@ -4766,7 -4791,7 +4766,7 @@@ dephii=0.0d0 dephii1=0.0d0 theti2=0.5d0*theta(i) - ityp2=ithetyp(iabs(itype(i-1))) + ityp2=ithetyp(itype(i-1)) do k=1,nntheterm coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) @@@ -4778,7 -4803,7 +4778,7 @@@ #else phii=phi(i) #endif - ityp1=ithetyp(iabs(itype(i-2))) + ityp1=ithetyp(itype(i-2)) do k=1,nsingle cosph1(k)=dcos(k*phii) sinph1(k)=dsin(k*phii) @@@ -4799,7 -4824,7 +4799,7 @@@ #else phii1=phi(i+1) #endif - ityp3=ithetyp(iabs(itype(i))) + ityp3=ithetyp(itype(i)) do k=1,nsingle cosph2(k)=dcos(k*phii1) sinph2(k)=dsin(k*phii1) @@@ -4950,7 -4975,7 +4950,7 @@@ c write (iout,'(a)') 'ESC do i=loc_start,loc_end it=itype(i) if (it.eq.10) goto 1 - nlobit=nlob(iabs(it)) + nlobit=nlob(it) c print *,'i=',i,' it=',it,' nlobit=',nlobit c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad theti=theta(i+1)-pipol @@@ -5107,11 -5132,11 +5107,11 @@@ C Compute the contribution to SC energ do j=1,nlobit #ifdef OSF - adexp=bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin + adexp=bsc(j,it)-0.5D0*contr(j,iii)+emin if(adexp.ne.adexp) adexp=1.0 expfac=dexp(adexp) #else - expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin) + expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin) #endif cd print *,'j=',j,' expfac=',expfac escloc_i=escloc_i+expfac @@@ -5193,7 -5218,7 +5193,7 @@@ C Compute the contribution to SC energ dersc12=0.0d0 do j=1,nlobit - expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j)+emin) + expfac=dexp(bsc(j,it)-0.5D0*contr(j)+emin) escloc_i=escloc_i+expfac do k=1,2 dersc(k)=dersc(k)+Ax(k,j)*expfac @@@ -5295,7 -5320,7 +5295,7 @@@ do j = 1,3 xx = xx + x_prime(j)*dc_norm(j,i+nres) yy = yy + y_prime(j)*dc_norm(j,i+nres) - zz = zz + z_prime(j)*dc_norm(j,i+nres) + zz = zz + dsign(1.0,itype(i))*z_prime(j)*dc_norm(j,i+nres) enddo xxtab(i)=xx @@@ -5305,7 -5330,7 +5305,7 @@@ C Compute the energy of the ith side cbain C c write (2,*) "xx",xx," yy",yy," zz",zz - it=itype(i) + it=iabs(itype(i)) do j = 1,65 x(j) = sc_parmin(j,it) enddo @@@ -5313,7 -5338,7 +5313,7 @@@ Cc diagnostics - remove later xx1 = dcos(alph(2)) yy1 = dsin(alph(2))*dcos(omeg(2)) - zz1 = -dsin(alph(2))*dsin(omeg(2)) + zz1 = -dsign(1.0,itype(i))*dsin(alph(2))*dsin(omeg(2)) write(2,'(3f8.1,3f9.3,1x,3f9.3)') & alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz, & xx1,yy1,zz1 @@@ -5772,12 -5797,17 +5772,12 @@@ c lprn=.true etors_ii=0.0D0 itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) phii=phi(i) gloci=0.0D0 C Regular cosine and sine terms - do j=1,nterm(itori,itori1,iblock) - v1ij=v1(j,itori,itori1,iblock) - v2ij=v2(j,itori,itori1,iblock) + do j=1,nterm(itori,itori1) + v1ij=v1(j,itori,itori1) + v2ij=v2(j,itori,itori1) cosphi=dcos(j*phii) sinphi=dsin(j*phii) etors=etors+v1ij*cosphi+v2ij*sinphi @@@ -5792,7 -5822,7 +5792,7 @@@ C [v2 cos(phi/2)+v3 sin(phi/2) C cosphi=dcos(0.5d0*phii) sinphi=dsin(0.5d0*phii) - do j=1,nlor(itori,itori1,iblock) + do j=1,nlor(itori,itori1) vl1ij=vlor1(j,itori,itori1) vl2ij=vlor2(j,itori,itori1) vl3ij=vlor3(j,itori,itori1) @@@ -5805,13 -5835,14 +5805,13 @@@ gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom enddo C Subtract the constant term - etors=etors-v0(itori,itori1,iblock) + etors=etors-v0(itori,itori1) if (energy_dec) write (iout,'(a6,i5,0pf7.3)') - & 'etor',i,etors_ii-v0(itori,itori1,iblock) + & 'etor',i,etors_ii-v0(itori,itori1) if (lprn) & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1, - & (v1(j,itori,itori1,iblock),j=1,6), - & (v2(j,itori,itori1,iblock),j=1,6) + & (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6) gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) enddo @@@ -5866,15 -5897,17 +5866,15 @@@ c lprn=.true itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) - iblock=1 - if (iabs(itype(i+1)).eq.20) iblock=2 phii=phi(i) phii1=phi(i+1) gloci1=0.0D0 gloci2=0.0D0 - do j=1,ntermd_1(itori,itori1,itori2,iblock) - v1cij=v1c(1,j,itori,itori1,itori2,iblock) - v1sij=v1s(1,j,itori,itori1,itori2,iblock) - v2cij=v1c(2,j,itori,itori1,itori2,iblock) - v2sij=v1s(2,j,itori,itori1,itori2,iblock) + do j=1,ntermd_1(itori,itori1,itori2) + v1cij=v1c(1,j,itori,itori1,itori2) + v1sij=v1s(1,j,itori,itori1,itori2) + v2cij=v1c(2,j,itori,itori1,itori2) + v2sij=v1s(2,j,itori,itori1,itori2) cosphi1=dcos(j*phii) sinphi1=dsin(j*phii) cosphi2=dcos(j*phii1) @@@ -5884,12 -5917,12 +5884,12 @@@ gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1) gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2) enddo - do k=2,ntermd_2(itori,itori1,itori2,iblock) + do k=2,ntermd_2(itori,itori1,itori2) do l=1,k-1 - v1cdij = v2c(k,l,itori,itori1,itori2,iblock) - v2cdij = v2c(l,k,itori,itori1,itori2,iblock) - v1sdij = v2s(k,l,itori,itori1,itori2,iblock) - v2sdij = v2s(l,k,itori,itori1,itori2,iblock) + v1cdij = v2c(k,l,itori,itori1,itori2) + v2cdij = v2c(l,k,itori,itori1,itori2) + v1sdij = v2s(k,l,itori,itori1,itori2) + v2sdij = v2s(l,k,itori,itori1,itori2) cosphi1p2=dcos(l*phii+(k-l)*phii1) cosphi1m2=dcos(l*phii-(k-l)*phii1) sinphi1p2=dsin(l*phii+(k-l)*phii1) @@@ -5939,6 -5972,7 +5939,7 @@@ c write (iout,*) "EBACK_SC_COR",ip esccor=0.0D0 do i=itau_start,itau_end esccor_ii=0.0D0 + if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle isccori=isccortyp(itype(i-2)) isccori1=isccortyp(itype(i-1)) phii=phi(i) @@@ -5957,14 -5991,14 +5958,14 @@@ c 2 = Ca...Ca...Ca...S c 3 = SC...Ca...Ca...SCi gloci=0.0D0 if (((intertyp.eq.3).and.((itype(i-2).eq.10).or. - & (itype(i-1).eq.10).or.(itype(i-2).eq.21).or. - & (itype(i-1).eq.21))) + & (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or. + & (itype(i-1).eq.ntyp1))) & .or. ((intertyp.eq.1).and.((itype(i-2).eq.10) - & .or.(itype(i-2).eq.21))) + & .or.(itype(i-2).eq.ntyp1))) & .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or. - & (itype(i-1).eq.21)))) cycle - if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.21)) cycle - if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.21)) + & (itype(i-1).eq.ntyp1)))) cycle + if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle + if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1)) & cycle do j=1,nterm_sccor(isccori,isccori1) v1ij=v1sccor(j,intertyp,isccori,isccori1) diff --combined source/unres/src_MD/parmread.F index 4729ca5,e15c7ae..54918a9 --- a/source/unres/src_MD/parmread.F +++ b/source/unres/src_MD/parmread.F @@@ -28,6 -28,7 +28,6 @@@ include 'COMMON.SETUP' character*1 t1,t2,t3 character*1 onelett(4) /"G","A","P","D"/ - character*1 toronelet(-2:2) /"p","a","G","A","P"/ logical lprint,LaTeX dimension blower(3,3,maxlob) dimension b(13) @@@ -102,13 -103,47 +102,13 @@@ C Read the parameters of the probabilit C of the virtual-bond valence angles theta C do i=1,ntyp - read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2), - & (bthet(j,i,1,1),j=1,2) + read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2), + & (bthet(j,i),j=1,2) read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3) read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3) read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i) sigc0(i)=sigc0(i)**2 enddo - do i=1,ntyp - athet(1,i,1,-1)=athet(1,i,1,1) - athet(2,i,1,-1)=athet(2,i,1,1) - bthet(1,i,1,-1)=-bthet(1,i,1,1) - bthet(2,i,1,-1)=-bthet(2,i,1,1) - athet(1,i,-1,1)=-athet(1,i,1,1) - athet(2,i,-1,1)=-athet(2,i,1,1) - bthet(1,i,-1,1)=bthet(1,i,1,1) - bthet(2,i,-1,1)=bthet(2,i,1,1) - enddo - do i=-ntyp,-1 - a0thet(i)=a0thet(-i) - athet(1,i,-1,-1)=athet(1,-i,1,1) - athet(2,i,-1,-1)=-athet(2,-i,1,1) - bthet(1,i,-1,-1)=bthet(1,-i,1,1) - bthet(2,i,-1,-1)=-bthet(2,-i,1,1) - athet(1,i,-1,1)=athet(1,-i,1,1) - athet(2,i,-1,1)=-athet(2,-i,1,1) - bthet(1,i,-1,1)=-bthet(1,-i,1,1) - bthet(2,i,-1,1)=bthet(2,-i,1,1) - athet(1,i,1,-1)=-athet(1,-i,1,1) - athet(2,i,1,-1)=athet(2,-i,1,1) - bthet(1,i,1,-1)=bthet(1,-i,1,1) - bthet(2,i,1,-1)=-bthet(2,-i,1,1) - theta0(i)=theta0(-i) - sig0(i)=sig0(-i) - sigc0(i)=sigc0(-i) - do j=0,3 - polthet(j,i)=polthet(j,-i) - enddo - do j=1,3 - gthet(j,i)=gthet(j,-i) - enddo - enddo close (ithep) if (lprint) then if (.not.LaTeX) then @@@ -119,7 -154,7 +119,7 @@@ & ' B1 ',' B2 ' do i=1,ntyp write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i, - & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2) + & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2) enddo write (iout,'(/a/9x,5a/79(1h-))') & 'Parameters of the expression for sigma(theta_c):', @@@ -146,7 -181,8 +146,7 @@@ & ' b1*10^1 ',' b2*10^1 ' do i=1,ntyp write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i), - & a0thet(i),(100*athet(j,i,1,1),j=1,2), - & (10*bthet(j,i,1,1),j=1,2) + & a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2) enddo write (iout,'(/a/9x,5a/79(1h-))') & 'Parameters of the expression for sigma(theta_c):', @@@ -312,10 -348,18 +312,10 @@@ bsc(1,i)=0.0D0 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3), & ((blower(k,l,1),l=1,k),k=1,3) - censc(1,1,-i)=censc(1,1,i) - censc(2,1,-i)=censc(2,1,i) - censc(3,1,-i)=-censc(3,1,i) - do j=2,nlob(i) read (irotam,*,end=112,err=112) bsc(j,i) read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3), & ((blower(k,l,j),l=1,k),k=1,3) - censc(1,j,-i)=censc(1,j,i) - censc(2,j,-i)=censc(2,j,i) - censc(3,j,-i)=-censc(3,j,i) -C BSC is amplitude of Gaussian enddo do j=1,nlob(i) do k=1,3 @@@ -326,6 -370,14 +326,6 @@@ enddo gaussc(k,l,j,i)=akl gaussc(l,k,j,i)=akl - if (((k.eq.3).and.(l.ne.3)) - & .or.((l.eq.3).and.(k.ne.3))) then - gaussc(k,l,j,-i)=-akl - gaussc(l,k,j,-i)=-akl - else - gaussc(k,l,j,-i)=akl - gaussc(l,k,j,-i)=akl - endif enddo enddo enddo @@@ -450,25 -502,40 +450,25 @@@ C Read torsional parameter C read (itorp,*,end=113,err=113) ntortyp read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp) - do iblock=1,2 - do i=-ntyp,-1 - itortyp(i)=-itortyp(-i) - enddo c write (iout,*) 'ntortyp',ntortyp - do i=0,ntortyp-1 - do j=-ntortyp+1,ntortyp-1 - read (itorp,*,end=113,err=113) nterm(i,j,iblock), - & nlor(i,j,iblock) - nterm(-i,-j,iblock)=nterm(i,j,iblock) - nlor(-i,-j,iblock)=nlor(i,j,iblock) + do i=1,ntortyp + do j=1,ntortyp + read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j) v0ij=0.0d0 si=-1.0d0 - do k=1,nterm(i,j,iblock) - read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock), - & v2(k,i,j,iblock) - v1(k,-i,-j,iblock)=v1(k,i,j,iblock) - v2(k,-i,-j,iblock)=-v2(k,i,j,iblock) - v0ij=v0ij+si*v1(k,i,j,iblock) + do k=1,nterm(i,j) + read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j) + v0ij=v0ij+si*v1(k,i,j) si=-si -c write(iout,*) i,j,k,iblock,nterm(i,j,iblock) -c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock), -c &v2(k,-i,-j,iblock),v2(k,i,j,iblock) enddo - do k=1,nlor(i,j,iblock) + do k=1,nlor(i,j) read (itorp,*,end=113,err=113) kk,vlor1(k,i,j), & vlor2(k,i,j),vlor3(k,i,j) v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2) enddo - v0(i,j,iblock)=v0ij - v0(-i,-j,iblock)=v0ij + v0(i,j)=v0ij enddo enddo - enddo close (itorp) if (lprint) then write (iout,'(/a/)') 'Torsional constants:' @@@ -476,11 -543,12 +476,11 @@@ do j=1,ntortyp write (iout,*) 'ityp',i,' jtyp',j write (iout,*) 'Fourier constants' - do k=1,nterm(i,j,iblock) - write (iout,'(2(1pe15.5))') v1(k,i,j,iblock), - & v2(k,i,j,iblock) + do k=1,nterm(i,j) + write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j) enddo write (iout,*) 'Lorenz constants' - do k=1,nlor(i,j,iblock) + do k=1,nlor(i,j) write (iout,'(3(1pe15.5))') & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j) enddo @@@ -490,12 -558,16 +490,15 @@@ C C 6/23/01 Read parameters for double torsionals C - do iblock=1,2 - do i=0,ntortyp-1 - do j=-ntortyp+1,ntortyp-1 - do k=-ntortyp+1,ntortyp-1 + do i=1,ntortyp + do j=1,ntortyp + do k=1,ntortyp read (itordp,'(3a1)',end=114,err=114) t1,t2,t3 - if (t1.ne.onelett(i) .or. t2.ne.onelett(j) - & .or. t3.ne.onelett(k)) then + c write (iout,*) "OK onelett", + c & i,j,k,t1,t2,t3 + + if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) + & .or. t3.ne.toronelet(k)) then write (iout,*) "Error in double torsional parameter file", & i,j,k,t1,t2,t3 #ifdef MPI @@@ -503,54 -575,81 +506,54 @@@ #endif stop "Error in double torsional parameter file" endif - read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock), - & ntermd_2(i,j,k,iblock) - ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock) - ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock) - read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1, - & ntermd_1(i,j,k,iblock)) - read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1, - & ntermd_1(i,j,k,iblock)) - read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1, - & ntermd_1(i,j,k,iblock)) - read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1, - & ntermd_1(i,j,k,iblock)) -C Martix of D parameters for one dimesional foureir series - do l=1,ntermd_1(i,j,k,iblock) - v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock) - v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock) - v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock) - v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock) -c write(iout,*) "whcodze" , -c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock) - enddo - read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock), - & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock), - & v2s(m,l,i,j,k,iblock), - & m=1,l-1),l=1,ntermd_2(i,j,k,iblock)) -C Martix of D parameters for two dimesional fourier series - do l=1,ntermd_2(i,j,k,iblock) - do m=1,l-1 - v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock) - v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock) - v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock) - v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock) - enddo!m - enddo!l - enddo!k - enddo!j - enddo!i - enddo!iblock + read (itordp,*,end=114,err=114) ntermd_1(i,j,k), + & ntermd_2(i,j,k) + read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1, + & ntermd_1(i,j,k)) + read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1, + & ntermd_1(i,j,k)) + read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1, + & ntermd_1(i,j,k)) + read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1, + & ntermd_1(i,j,k)) + read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k), + & v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k), + & m=1,l-1),l=1,ntermd_2(i,j,k)) + enddo + enddo + enddo if (lprint) then write (iout,*) write (iout,*) 'Constants for double torsionals' - do iblock=1,2 do i=1,ntortyp - do j=-ntortyp,ntortyp - do k=-ntortyp,ntortyp + do j=1,ntortyp + do k=1,ntortyp write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k, - & ' nsingle',ntermd_1(i,j,k,iblock), - & ' ndouble',ntermd_2(i,j,k,iblock) + & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k) write (iout,*) write (iout,*) 'Single angles:' - do l=1,ntermd_1(i,j,k,iblock) - write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l, - & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock), - & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock), - & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock) + do l=1,ntermd_1(i,j,k) + write (iout,'(i5,2f10.5,5x,2f10.5)') l, + & v1c(1,l,i,j,k),v1s(1,l,i,j,k), + & v1c(2,l,i,j,k),v1s(2,l,i,j,k) enddo write (iout,*) write (iout,*) 'Pairs of angles:' - write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock)) - do l=1,ntermd_2(i,j,k,iblock) + write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k)) + do l=1,ntermd_2(i,j,k) write (iout,'(i5,20f10.5)') - & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)) + & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k)) enddo write (iout,*) - write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock)) - do l=1,ntermd_2(i,j,k,iblock) + write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k)) + do l=1,ntermd_2(i,j,k) write (iout,'(i5,20f10.5)') - & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)), - & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock)) + & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k)) enddo write (iout,*) enddo enddo enddo - enddo endif #endif C Read of Side-chain backbone correlation parameters @@@ -558,7 -657,12 +561,11 @@@ C Modified 11 May 2012 by Adask CCC C read (isccor,*,end=113,err=113) nsccortyp -#ifdef SCCORPDB read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp) + do i=-ntyp,-1 + isccortyp(i)=-isccortyp(-i) + enddo + iscprol=isccortyp(20) c write (iout,*) 'ntortyp',ntortyp maxinter=3 cc maxinter is maximum interaction sites @@@ -567,12 -671,52 +574,52 @@@ do j=1,nsccortyp read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j) v0ijsccor=0.0d0 + v0ijsccor1=0.0d0 + v0ijsccor2=0.0d0 + v0ijsccor3=0.0d0 si=-1.0d0 - + nterm_sccor(-i,j)=nterm_sccor(i,j) + nterm_sccor(-i,-j)=nterm_sccor(i,j) + nterm_sccor(i,-j)=nterm_sccor(i,j) do k=1,nterm_sccor(i,j) read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j) - & ,v2sccor(k,l,i,j) + & ,v2sccor(k,l,i,j) + if (j.eq.iscprol) then + if (i.eq.isccortyp(10)) then + v1sccor(k,l,i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) + else + v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 + & +v2sccor(k,l,i,j)*dsqrt(0.75d0) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 + & +v1sccor(k,l,i,j)*dsqrt(0.75d0) + v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j) + v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j) + v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j) + endif + else + if (i.eq.isccortyp(10)) then + v1sccor(k,l,i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) + else + if (j.eq.isccortyp(10)) then + v1sccor(k,l,-i,j)=v1sccor(k,l,i,j) + v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j) + else + v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) + v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j) + v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j) + v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j) + endif + endif + endif v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j) + v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j) + v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j) + v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j) si=-si enddo do k=1,nlor_sccor(i,j) @@@ -581,12 -725,46 +628,46 @@@ v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &(1+vlor3sccor(k,i,j)**2) enddo + v0sccor(l,i,j)=v0ijsccor + v0sccor(l,-i,j)=v0ijsccor1 + v0sccor(l,i,-j)=v0ijsccor2 + v0sccor(l,-i,-j)=v0ijsccor3 + enddo + enddo + enddo + close (isccor) + #else + read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp) + c write (iout,*) 'ntortyp',ntortyp + maxinter=3 + cc maxinter is maximum interaction sites + do l=1,maxinter + do i=1,nsccortyp + do j=1,nsccortyp + read (isccor,*,end=113,err=113) + & nterm_sccor(i,j),nlor_sccor(i,j) + v0ijsccor=0.0d0 + si=-1.0d0 + + do k=1,nterm_sccor(i,j) + read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j) + & ,v2sccor(k,l,i,j) + v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j) + si=-si + enddo + do k=1,nlor_sccor(i,j) + read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j), + & vlor2sccor(k,i,j),vlor3sccor(k,i,j) + v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ + &(1+vlor3sccor(k,i,j)**2) + enddo v0sccor(i,j)=v0ijsccor enddo enddo enddo close (isccor) - + + #endif if (lprint) then write (iout,'(/a/)') 'Torsional constants:' do i=1,nsccortyp @@@ -614,7 -792,7 +695,7 @@@ write (iout,*) "Coefficients of the cumulants" endif read (ifourier,*) nloctyp - do i=0,nloctyp-1 + do i=1,nloctyp read (ifourier,*,end=115,err=115) read (ifourier,*,end=115,err=115) (b(ii),ii=1,13) if (lprint) then @@@ -623,20 -801,31 +704,22 @@@ endif B1(1,i) = b(3) B1(2,i) = b(5) c b1(1,i)=0.0d0 c b1(2,i)=0.0d0 B1tilde(1,i) = b(3) - B1tilde(2,i) =-b(5) + B1tilde(2,i) =-b(5) + B1tilde(1,-i) =-b(3) + B1tilde(2,-i) =b(5) c b1tilde(1,i)=0.0d0 c b1tilde(2,i)=0.0d0 B2(1,i) = b(2) B2(2,i) = b(4) - B2(1,-i) =b(2) - B2(2,-i) =-b(4) - c b2(1,i)=0.0d0 c b2(2,i)=0.0d0 CC(1,1,i)= b(7) CC(2,2,i)=-b(7) CC(2,1,i)= b(9) CC(1,2,i)= b(9) - CC(1,1,-i)= b(7) - CC(2,2,-i)=-b(7) - CC(2,1,-i)=-b(9) - CC(1,2,-i)=-b(9) c CC(1,1,i)=0.0d0 c CC(2,2,i)=0.0d0 c CC(2,1,i)=0.0d0 @@@ -645,6 -834,11 +728,6 @@@ c CC(1,2,i)=0.0d Ctilde(1,2,i)=b(9) Ctilde(2,1,i)=-b(9) Ctilde(2,2,i)=b(7) - Ctilde(1,1,-i)=b(7) - Ctilde(1,2,-i)=-b(9) - Ctilde(2,1,-i)=b(9) - Ctilde(2,2,-i)=b(7) - c Ctilde(1,1,i)=0.0d0 c Ctilde(1,2,i)=0.0d0 c Ctilde(2,1,i)=0.0d0 @@@ -653,6 -847,10 +736,6 @@@ c Ctilde(2,2,i)=0.0d DD(2,2,i)=-b(6) DD(2,1,i)= b(8) DD(1,2,i)= b(8) - DD(1,1,-i)= b(6) - DD(2,2,-i)=-b(6) - DD(2,1,-i)=-b(8) - DD(1,2,-i)=-b(8) c DD(1,1,i)=0.0d0 c DD(2,2,i)=0.0d0 c DD(2,1,i)=0.0d0 @@@ -661,6 -859,11 +744,6 @@@ c DD(1,2,i)=0.0d Dtilde(1,2,i)=b(8) Dtilde(2,1,i)=-b(8) Dtilde(2,2,i)=b(6) - Dtilde(1,1,-i)=b(6) - Dtilde(1,2,-i)=-b(8) - Dtilde(2,1,-i)=b(8) - Dtilde(2,2,-i)=b(6) - c Dtilde(1,1,i)=0.0d0 c Dtilde(1,2,i)=0.0d0 c Dtilde(2,1,i)=0.0d0 @@@ -669,6 -872,11 +752,6 @@@ c Dtilde(2,2,i)=0.0d EE(2,2,i)=-b(10)+b(11) EE(2,1,i)= b(12)-b(13) EE(1,2,i)= b(12)+b(13) - EE(1,1,-i)= b(10)+b(11) - EE(2,2,-i)=-b(10)+b(11) - EE(2,1,-i)=-b(12)+b(13) - EE(1,2,-i)=-b(12)-b(13) - c ee(1,1,i)=1.0d0 c ee(2,2,i)=1.0d0 c ee(2,1,i)=0.0d0 diff --combined source/unres/src_MD/readrtns.F index 47850c2,25a4bc0..c15396a --- a/source/unres/src_MD/readrtns.F +++ b/source/unres/src_MD/readrtns.F @@@ -930,8 -930,8 +930,8 @@@ C Assign initial virtual bond length vbld_inv(i)=vblinv enddo do i=2,nres-1 - vbld(i+nres)=dsc(iabs(itype(i))) - vbld_inv(i+nres)=dsc_inv(iabs(itype(i))) + vbld(i+nres)=dsc(itype(i)) + vbld_inv(i+nres)=dsc_inv(itype(i)) c write (iout,*) "i",i," itype",itype(i), c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres) enddo @@@ -940,15 -940,15 +940,15 @@@ c print *,nre c print '(20i4)',(itype(i),i=1,nres) do i=1,nres #ifdef PROCOR - if (itype(i).eq.21 .or. itype(i+1).eq.21) then + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then #else - if (itype(i).eq.21) then + if (itype(i).eq.ntyp1) then #endif itel(i)=0 #ifdef PROCOR - else if (iabs(itype(i+1)).ne.20) then + else if (itype(i+1).ne.20) then #else - else if (iabs(itype(i)).ne.20) then + else if (itype(i).ne.20) then #endif itel(i)=1 else @@@ -1005,8 -1005,8 +1005,8 @@@ C 8/13/98 Set limits to generating the #endif nct=nres cd print *,'NNT=',NNT,' NCT=',NCT - if (itype(1).eq.21) nnt=2 - if (itype(nres).eq.21) nct=nct-1 + if (itype(1).eq.ntyp1) nnt=2 + if (itype(nres).eq.ntyp1) nct=nct-1 if (pdbref) then if(me.eq.king.or..not.out1file) & write (iout,'(a,i3)') 'nsup=',nsup @@@ -1151,6 -1151,7 +1151,6 @@@ C initial geometry enddo do i=2,nres-1 omeg(i)=-120d0*deg2rad - if (itype(i).le.0) omeg(i)=-omeg(i) enddo else if(me.eq.king.or..not.out1file) diff --combined source/unres/src_MD/sc_move.F index b6837fd,274767b..02a329d --- a/source/unres/src_MD/sc_move.F +++ b/source/unres/src_MD/sc_move.F @@@ -213,7 -213,7 +213,7 @@@ c Define what is meant by "neighbou c Don't do glycine or ends i=itype(res_pick) - if (i.eq.10 .or. i.eq.21) return + if (i.eq.10 .or. i.eq.ntyp1) return c Freeze everything (later will relax only selected side-chains) mask_r=.true. @@@ -255,7 -255,7 +255,7 @@@ cd print *,'new ',(energy(k) n_try=0 do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop) c Move the selected residue (don't worry if it fails) - call gen_side(iabs(itype(res_pick)),theta(res_pick+1), + call gen_side(itype(res_pick),theta(res_pick+1), + alph(res_pick),omeg(res_pick),fail) c Minimize the side-chains starting from the new arrangement @@@ -719,8 -719,8 +719,8 @@@ c if (icall.eq.0) lprn=.true do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) - itypi1=iabs(itype(i+1)) + itypi=itype(i) + itypi1=itype(i+1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@@ -735,7 -735,7 +735,7 @@@ do j=istart(i,iint),iend(i,iint) IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN ind=ind+1 - itypj=iabs(itype(j)) + itypj=itype(j) dscj_inv=dsc_inv(itypj) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) diff --combined source/wham/src/enecalc1.F index c431153,ff45881..3695ba7 --- a/source/wham/src/enecalc1.F +++ b/source/wham/src/enecalc1.F @@@ -195,6 -195,11 +195,10 @@@ c call pdbout(ii+1,beta_h(ib,ipa & iii+1,indstart(me1)+iii," T", & 1.0d0/(1.987D-3*beta_h(ib,ipar)) errmsg_count=errmsg_count+1 - + call pdbout(indstart(me1)+iii, + & 1.0d0/(1.987D-3*beta_h(ib,ipar)), + &energia(0),eini,0.0d0,0.0d0) + call enerprint(energia(0),fT) if (errmsg_count.gt.maxerrmsg_count) & write (iout,*) "Too many warning messages" if (einicheck.gt.1) then @@@ -725,7 -730,8 +729,8 @@@ c-------------------------------------- enddo do j=nnt,nct itj=itype(j) - if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(itj)).gt.2.0d0) then + if (itype(j).ne.10 .and.(vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) + & then if (iprint.gt.0) & write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j), & " for conformation",ii diff --combined source/wham/src/energy_p_new.F index cb4b7b0,8303bc2..2b48cf6 --- a/source/wham/src/energy_p_new.F +++ b/source/wham/src/energy_p_new.F @@@ -368,8 -368,8 +368,8 @@@ cd print *,'Entering ELJ nnt=',nnt, evdw=0.0D0 evdw_t=0.0d0 do i=iatsc_s,iatsc_e - itypi=itype(i) - itypi1=itype(i+1) + itypi=iabs(itype(i)) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@@ -382,7 -382,7 +382,7 @@@ cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) + itypj=iabs(itype(j)) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@@ -539,8 -539,8 +539,8 @@@ c print *,'Entering ELJK nnt=',nnt, evdw=0.0D0 evdw_t=0.0d0 do i=iatsc_s,iatsc_e - itypi=itype(i) - itypi1=itype(i+1) + itypi=iabs(itype(i)) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@@ -549,7 -549,7 +549,7 @@@ C Calculate SC interaction energy C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) + itypj=iabs(itype(j)) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi zj=c(3,nres+j)-zi @@@ -650,8 -650,8 +650,8 @@@ c els c endif ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) - itypi1=itype(i+1) + itypi=iabs(itype(i)) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@@ -665,7 -665,7 +665,7 @@@ do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) dscj_inv=vbld_inv(j+nres) chi1=chi(itypi,itypj) chi2=chi(itypj,itypi) @@@ -786,8 -786,8 +786,8 @@@ c print *,'Entering EGB nnt=',nnt, c if (icall.gt.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) - itypi1=itype(i+1) + itypi=iabs(itype(i)) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@@ -801,7 -801,7 +801,7 @@@ do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) @@@ -931,8 -931,8 +931,8 @@@ c print *,'Entering EGB nnt=',nnt, c if (icall.gt.0) lprn=.true. ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) - itypi1=itype(i+1) + itypi=iabs(itype(i)) + itypi1=iabs(itype(i+1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@@ -946,7 -946,7 +946,7 @@@ do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) ind=ind+1 - itypj=itype(j) + itypj=iabs(itype(j)) dscj_inv=vbld_inv(j+nres) sig0ij=sigma(itypi,itypj) r0ij=r0(itypi,itypj) @@@ -2785,7 -2785,7 +2785,7 @@@ c & " iscp",(iscpstart(i,j),iscpe do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) - itypj=itype(j) + itypj=iabs(itype(j)) C Uncomment following three lines for SC-p interactions c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi @@@ -2897,7 -2897,8 +2897,8 @@@ c write (iout,*) "i",i," ii",ii, c & dhpb(i),dhpb1(i),forcon(i) C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. - if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then + if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. + & iabs(itype(jjj)).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij cd write (iout,*) "eij",eij @@@ -2997,7 -2998,7 +2998,7 @@@ include 'COMMON.VAR' include 'COMMON.IOUNITS' double precision erij(3),dcosom1(3),dcosom2(3),gg(3) - itypi=itype(i) + itypi=iabs(itype(i)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@@ -3005,7 -3006,7 +3006,7 @@@ dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) dsci_inv=dsc_inv(itypi) - itypj=itype(j) + itypj=iabs(itype(j)) dscj_inv=dsc_inv(itypj) xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi @@@ -3095,7 -3096,7 +3096,7 @@@ c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included c do i=nnt,nct - iti=itype(i) + iti=iabs(itype(i)) if (iti.ne.10) then nbi=nbondterm(iti) if (nbi.eq.1) then @@@ -3175,6 -3176,18 +3176,18 @@@ c write (iout,*) ithet_start,ithet C Zero the energy function and its derivative at 0 or pi. call splinthet(theta(i),0.5d0*delta,ss,ssd) it=itype(i-1) + ichir1=isign(1,itype(i-2)) + ichir2=isign(1,itype(i)) + if (itype(i-2).eq.10) ichir1=isign(1,itype(i-1)) + if (itype(i).eq.10) ichir2=isign(1,itype(i-1)) + if (itype(i-1).eq.10) then + itype1=isign(10,itype(i-2)) + ichir11=isign(1,itype(i-2)) + ichir12=isign(1,itype(i-2)) + itype2=isign(10,itype(i)) + ichir21=isign(1,itype(i)) + ichir22=isign(1,itype(i)) + endif c if (i.gt.ithet_start .and. c & (itel(i-1).eq.0 .or. itel(i-2).eq.0)) goto 1215 c if (i.gt.3 .and. (i.le.4 .or. itel(i-3).ne.0)) then @@@ -3230,8 -3243,12 +3243,12 @@@ C dependent on the adjacent virtual-bon C In following comments this theta will be referred to as t_c. thet_pred_mean=0.0d0 do k=1,2 - athetk=athet(k,it) - bthetk=bthet(k,it) + athetk=athet(k,it,ichir1,ichir2) + bthetk=bthet(k,it,ichir1,ichir2) + if (it.eq.10) then + athetk=athet(k,itype1,ichir11,ichir12) + bthetk=bthet(k,itype2,ichir21,ichir22) + endif thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k) enddo c write (iout,*) "thet_pred_mean",thet_pred_mean @@@ -3239,8 -3256,16 +3256,16 @@@ thet_pred_mean=thet_pred_mean*ss+a0thet(it) c write (iout,*) "thet_pred_mean",thet_pred_mean C Derivatives of the "mean" values in gamma1 and gamma2. - dthetg1=(-athet(1,it)*y(2)+athet(2,it)*y(1))*ss - dthetg2=(-bthet(1,it)*z(2)+bthet(2,it)*z(1))*ss + dthetg1=(-athet(1,it,ichir1,ichir2)*y(2) + &+athet(2,it,ichir1,ichir2)*y(1))*ss + dthetg2=(-bthet(1,it,ichir1,ichir2)*z(2) + & +bthet(2,it,ichir1,ichir2)*z(1))*ss + if (it.eq.10) then + dthetg1=(-athet(1,itype1,ichir11,ichir12)*y(2) + &+athet(2,itype1,ichir11,ichir12)*y(1))*ss + dthetg2=(-bthet(1,itype2,ichir21,ichir22)*z(2) + & +bthet(2,itype2,ichir21,ichir22)*z(1))*ss + endif if (theta(i).gt.pi-delta) then call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0, & E_tc0) @@@ -3415,7 -3440,7 +3440,7 @@@ c write (iout,*) "ithetyp",(ithety dephii=0.0d0 dephii1=0.0d0 theti2=0.5d0*theta(i) - ityp2=ithetyp(itype(i-1)) + ityp2=ithetyp(iabs(itype(i-1))) do k=1,nntheterm coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) @@@ -3427,7 -3452,7 +3452,7 @@@ #else phii=phi(i) #endif - ityp1=ithetyp(itype(i-2)) + ityp1=ithetyp(iabs(itype(i-2))) do k=1,nsingle cosph1(k)=dcos(k*phii) sinph1(k)=dsin(k*phii) @@@ -3448,7 -3473,7 +3473,7 @@@ #else phii1=phi(i+1) #endif - ityp3=ithetyp(itype(i)) + ityp3=ithetyp(iabs(itype(i))) do k=1,nsingle cosph2(k)=dcos(k*phii1) sinph2(k)=dsin(k*phii1) @@@ -3602,7 -3627,7 +3627,7 @@@ c write (iout,'(a)') 'ESC do i=loc_start,loc_end it=itype(i) if (it.eq.10) goto 1 - nlobit=nlob(it) + nlobit=nlob(iabs(it)) c print *,'i=',i,' it=',it,' nlobit=',nlobit c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad theti=theta(i+1)-pipol @@@ -3757,7 -3782,7 +3782,7 @@@ C Compute the contribution to SC energ do iii=-1,1 do j=1,nlobit - expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin) + expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin) cd print *,'j=',j,' expfac=',expfac escloc_i=escloc_i+expfac do k=1,3 @@@ -3838,7 -3863,7 +3863,7 @@@ C Compute the contribution to SC energ dersc12=0.0d0 do j=1,nlobit - expfac=dexp(bsc(j,it)-0.5D0*contr(j)+emin) + expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j)+emin) escloc_i=escloc_i+expfac do k=1,2 dersc(k)=dersc(k)+Ax(k,j)*expfac @@@ -3901,7 -3926,7 +3926,7 @@@ cosfac=dsqrt(cosfac2) sinfac2=0.5d0/(1.0d0-costtab(i+1)) sinfac=dsqrt(sinfac2) - it=itype(i) + it=iabs(itype(i)) if (it.eq.10) goto 1 c C Compute the axes of tghe local cartesian coordinates system; store in @@@ -3941,7 -3966,7 +3966,7 @@@ do j = 1,3 xx = xx + x_prime(j)*dc_norm(j,i+nres) yy = yy + y_prime(j)*dc_norm(j,i+nres) - zz = zz + z_prime(j)*dc_norm(j,i+nres) + zz = zz + dsign(1.0,itype(i))*z_prime(j)*dc_norm(j,i+nres) enddo xxtab(i)=xx @@@ -3951,7 -3976,7 +3976,7 @@@ C Compute the energy of the ith side cbain C c write (2,*) "xx",xx," yy",yy," zz",zz - it=itype(i) + it=iabs(itype(i)) do j = 1,65 x(j) = sc_parmin(j,it) enddo @@@ -3959,7 -3984,7 +3984,7 @@@ Cc diagnostics - remove later xx1 = dcos(alph(2)) yy1 = dsin(alph(2))*dcos(omeg(2)) - zz1 = -dsin(alph(2))*dsin(omeg(2)) + zz1 = -dsign(1.0,itype(i))*dsin(alph(2))*dsin(omeg(2)) write(2,'(3f8.1,3f9.3,1x,3f9.3)') & alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz, & xx1,yy1,zz1 @@@ -4367,14 -4392,19 +4392,19 @@@ c lprn=.true etors=0.0D0 do i=iphi_start,iphi_end if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215 + if (iabs(itype(i)).eq.20) then + iblock=2 + else + iblock=1 + endif itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) phii=phi(i) gloci=0.0D0 C Regular cosine and sine terms - do j=1,nterm(itori,itori1) - v1ij=v1(j,itori,itori1) - v2ij=v2(j,itori,itori1) + do j=1,nterm(itori,itori1,iblock) + v1ij=v1(j,itori,itori1,iblock) + v2ij=v2(j,itori,itori1,iblock) cosphi=dcos(j*phii) sinphi=dsin(j*phii) etors=etors+v1ij*cosphi+v2ij*sinphi @@@ -4387,7 -4417,7 +4417,7 @@@ C [v2 cos(phi/2)+v3 sin(phi/2) C cosphi=dcos(0.5d0*phii) sinphi=dsin(0.5d0*phii) - do j=1,nlor(itori,itori1) + do j=1,nlor(itori,itori1,iblock) vl1ij=vlor1(j,itori,itori1) vl2ij=vlor2(j,itori,itori1) vl3ij=vlor3(j,itori,itori1) @@@ -4398,11 -4428,11 +4428,11 @@@ gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom enddo C Subtract the constant term - etors=etors-v0(itori,itori1) + etors=etors-v0(itori,itori1,iblock) if (lprn) & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1, - & (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6) + & (v1(j,itori,itori1,1),j=1,6),(v2(j,itori,itori1,1),j=1,6) gloc(i-3,icg)=gloc(i-3,icg)+wtor*fact*gloci c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) 1215 continue @@@ -4463,16 -4493,26 +4493,23 @@@ c lprn=.true itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) -c iblock=1 -c if (iabs(itype(i+1)).eq.20) iblock=2 phii=phi(i) phii1=phi(i+1) gloci1=0.0D0 gloci2=0.0D0 + iblock=1 + if (iabs(itype(i+1)).eq.20) iblock=2 C Regular cosine and sine terms - do j=1,ntermd_1(itori,itori1,itori2) - v1cij=v1c(1,j,itori,itori1,itori2) - v1sij=v1s(1,j,itori,itori1,itori2) - v2cij=v1c(2,j,itori,itori1,itori2) - v2sij=v1s(2,j,itori,itori1,itori2) + c c do j=1,ntermd_1(itori,itori1,itori2,iblock) + c v1cij=v1c(1,j,itori,itori1,itori2,iblock) + c v1sij=v1s(1,j,itori,itori1,itori2,iblock) + c v2cij=v1c(2,j,itori,itori1,itori2,iblock) + c v2sij=v1s(2,j,itori,itori1,itori2,iblock) + do j=1,ntermd_1(itori,itori1,itori2,iblock) + v1cij=v1c(1,j,itori,itori1,itori2,iblock) + v1sij=v1s(1,j,itori,itori1,itori2,iblock) + v2cij=v1c(2,j,itori,itori1,itori2,iblock) + v2sij=v1s(2,j,itori,itori1,itori2,iblock) - cosphi1=dcos(j*phii) sinphi1=dsin(j*phii) cosphi2=dcos(j*phii1) @@@ -4482,12 -4522,12 +4519,12 @@@ gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1) gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2) enddo - do k=2,ntermd_2(itori,itori1,itori2) + do k=2,ntermd_2(itori,itori1,itori2,iblock) do l=1,k-1 - v1cdij = v2c(k,l,itori,itori1,itori2) - v2cdij = v2c(l,k,itori,itori1,itori2) - v1sdij = v2s(k,l,itori,itori1,itori2) - v2sdij = v2s(l,k,itori,itori1,itori2) + v1cdij = v2c(k,l,itori,itori1,itori2,iblock) + v2cdij = v2c(l,k,itori,itori1,itori2,iblock) + v1sdij = v2s(k,l,itori,itori1,itori2,iblock) + v2sdij = v2s(l,k,itori,itori1,itori2,iblock) cosphi1p2=dcos(l*phii+(k-l)*phii1) cosphi1m2=dcos(l*phii-(k-l)*phii1) sinphi1p2=dsin(l*phii+(k-l)*phii1) @@@ -4534,12 -4574,12 +4571,12 @@@ c amino-acid residues C Set lprn=.true. for debugging lprn=.false. c lprn=.true. -c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor +c write (iout,*) "EBACK_SC_COR",itau_start,itau_end,nterm_sccor esccor=0.0D0 do i=itau_start,itau_end esccor_ii=0.0D0 - isccori=isccortyp(itype(i-2)) - isccori1=isccortyp(itype(i-1)) + isccori=isccortyp((itype(i-2))) + isccori1=isccortyp((itype(i-1))) phii=phi(i) cccc Added 9 May 2012 cc Tauangle is torsional engle depending on the value of first digit @@@ -4556,14 -4596,14 +4593,14 @@@ c 2 = Ca...Ca...Ca...S c 3 = SC...Ca...Ca...SCi gloci=0.0D0 if (((intertyp.eq.3).and.((itype(i-2).eq.10).or. - & (itype(i-1).eq.10).or.(itype(i-2).eq.21).or. - & (itype(i-1).eq.21))) + & (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or. + & (itype(i-1).eq.ntyp1))) & .or. ((intertyp.eq.1).and.((itype(i-2).eq.10) - & .or.(itype(i-2).eq.21))) + & .or.(itype(i-2).eq.ntyp1))) & .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or. - & (itype(i-1).eq.21)))) cycle - if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.21)) cycle - if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.21)) + & (itype(i-1).eq.ntyp1)))) cycle + if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle + if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1)) & cycle do j=1,nterm_sccor(isccori,isccori1) v1ij=v1sccor(j,intertyp,isccori,isccori1) diff --combined source/wham/src/parmread.F index 1b66ce0,29f1021..1c1c818 --- a/source/wham/src/parmread.F +++ b/source/wham/src/parmread.F @@@ -23,6 -23,7 +23,6 @@@ include 'COMMON.FREE' character*1 t1,t2,t3 character*1 onelett(4) /"G","A","P","D"/ - character*1 toronelet(-2:2)/"p","a","G","A","P"/ logical lprint dimension blower(3,3,maxlob) character*800 controlcard @@@ -75,7 -76,6 +75,7 @@@ wtor=ww(13) wtor_d=ww(14) wvdwpp=ww(16) + wstrain=ww(15) wbond=ww(18) wsccor=ww(19) @@@ -194,12 -194,47 +194,47 @@@ C Read the parameters of the probabilit C of the virtual-bond valence angles theta C do i=1,ntyp - read (ithep,*) a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2) + read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2), + & (bthet(j,i,1,1),j=1,2) read (ithep,*) (polthet(j,i),j=0,3) read (ithep,*) (gthet(j,i),j=1,3) read (ithep,*) theta0(i),sig0(i),sigc0(i) sigc0(i)=sigc0(i)**2 enddo + do i=1,ntyp + athet(1,i,1,-1)=athet(1,i,1,1) + athet(2,i,1,-1)=athet(2,i,1,1) + bthet(1,i,1,-1)=-bthet(1,i,1,1) + bthet(2,i,1,-1)=-bthet(2,i,1,1) + athet(1,i,-1,1)=-athet(1,i,1,1) + athet(2,i,-1,1)=-athet(2,i,1,1) + bthet(1,i,-1,1)=bthet(1,i,1,1) + bthet(2,i,-1,1)=bthet(2,i,1,1) + enddo + do i=-ntyp,-1 + a0thet(i)=a0thet(-i) + athet(1,i,-1,-1)=athet(1,-i,1,1) + athet(2,i,-1,-1)=-athet(2,-i,1,1) + bthet(1,i,-1,-1)=bthet(1,-i,1,1) + bthet(2,i,-1,-1)=-bthet(2,-i,1,1) + athet(1,i,-1,1)=athet(1,-i,1,1) + athet(2,i,-1,1)=-athet(2,-i,1,1) + bthet(1,i,-1,1)=-bthet(1,-i,1,1) + bthet(2,i,-1,1)=bthet(2,-i,1,1) + athet(1,i,1,-1)=-athet(1,-i,1,1) + athet(2,i,1,-1)=athet(2,-i,1,1) + bthet(1,i,1,-1)=bthet(1,-i,1,1) + bthet(2,i,1,-1)=-bthet(2,-i,1,1) + theta0(i)=theta0(-i) + sig0(i)=sig0(-i) + sigc0(i)=sigc0(-i) + do j=0,3 + polthet(j,i)=polthet(j,-i) + enddo + do j=1,3 + gthet(j,i)=gthet(j,-i) + enddo + enddo close (ithep) if (lprint) then c write (iout,'(a)') @@@ -235,7 -270,8 +270,8 @@@ c endd & ' b1*10^1 ',' b2*10^1 ' do i=1,ntyp write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i), - & a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2) + & a0thet(i),(100*athet(j,i,1,1),j=1,2), + & (10*bthet(j,i,1,1),j=1,2) enddo write (iout,'(/a/9x,5a/79(1h-))') & 'Parameters of the expression for sigma(theta_c):', @@@ -388,10 -424,16 +424,16 @@@ enddo bsc(1,i)=0.0D0 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3) + censc(1,1,-i)=censc(1,1,i) + censc(2,1,-i)=censc(2,1,i) + censc(3,1,-i)=-censc(3,1,i) do j=2,nlob(i) read (irotam,*) bsc(j,i) read (irotam,*) (censc(k,j,i),k=1,3), & ((blower(k,l,j),l=1,k),k=1,3) + censc(1,j,-i)=censc(1,j,i) + censc(2,j,-i)=censc(2,j,i) + censc(3,j,-i)=-censc(3,j,i) enddo do j=1,nlob(i) do k=1,3 @@@ -402,6 -444,14 +444,14 @@@ enddo gaussc(k,l,j,i)=akl gaussc(l,k,j,i)=akl + if (((k.eq.3).and.(l.ne.3)) + & .or.((l.eq.3).and.(k.ne.3))) then + gaussc(k,l,j,-i)=-akl + gaussc(l,k,j,-i)=-akl + else + gaussc(k,l,j,-i)=akl + gaussc(l,k,j,-i)=akl + endif enddo enddo enddo @@@ -490,24 -540,35 +540,35 @@@ C Read torsional parameter C read (itorp,*) ntortyp read (itorp,*) (itortyp(i),i=1,ntyp) - write (iout,*) 'ntortyp',ntortyp - do i=1,ntortyp - do j=1,ntortyp - read (itorp,*) nterm(i,j),nlor(i,j) + do iblock=1,2 + do i=-ntyp,-1 + itortyp(i)=-itortyp(-i) + enddo + c write (iout,*) 'ntortyp',ntortyp + do i=0,ntortyp-1 + do j=-ntortyp+1,ntortyp-1 + read (itorp,*) nterm(i,j,iblock), + & nlor(i,j,iblock) + nterm(-i,-j,iblock)=nterm(i,j,iblock) + nlor(-i,-j,iblock)=nlor(i,j,iblock) v0ij=0.0d0 si=-1.0d0 - do k=1,nterm(i,j) - read (itorp,*) kk,v1(k,i,j),v2(k,i,j) - v0ij=v0ij+si*v1(k,i,j) + do k=1,nterm(i,j,iblock) + read (itorp,*) kk,v1(k,i,j,iblock),v2(k,i,j,iblock) + v1(k,-i,-j,iblock)=v1(k,i,j,iblock) + v2(k,-i,-j,iblock)=-v2(k,i,j,iblock) + v0ij=v0ij+si*v1(k,i,j,iblock) si=-si enddo - do k=1,nlor(i,j) + do k=1,nlor(i,j,iblock) read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j) v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2) enddo - v0(i,j)=v0ij + v0(i,j,iblock)=v0ij + v0(-i,-j,iblock)=v0ij enddo enddo + enddo close (itorp) if (lprint) then write (iout,'(/a/)') 'Torsional constants:' @@@ -515,11 -576,12 +576,12 @@@ do j=1,ntortyp write (iout,*) 'ityp',i,' jtyp',j write (iout,*) 'Fourier constants' - do k=1,nterm(i,j) - write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j) + do k=1,nterm(i,j,iblock) + write (iout,'(2(1pe15.5))') v1(k,i,j,iblock), + & v2(k,i,j,iblock) enddo write (iout,*) 'Lorenz constants' - do k=1,nlor(i,j) + do k=1,nlor(i,j,iblock) write (iout,'(3(1pe15.5))') & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j) enddo @@@ -529,58 -591,90 +591,90 @@@ C C 6/23/01 Read parameters for double torsionals C - do i=1,ntortyp - do j=1,ntortyp - do k=1,ntortyp + do iblock=1,2 + do i=0,ntortyp-1 + do j=-ntortyp+1,ntortyp-1 + do k=-ntortyp+1,ntortyp-1 read (itordp,'(3a1)') t1,t2,t3 - if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) - & .or. t3.ne.toronelet(k)) then + if (t1.ne.onelett(i) .or. t2.ne.onelett(j) + & .or. t3.ne.onelett(k)) then write (iout,*) "Error in double torsional parameter file", & i,j,k,t1,t2,t3 stop "Error in double torsional parameter file" endif - read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k) - read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k)) - read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k)) - read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k)) - read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k)) - read (itordp,*) ((v2c(l,m,i,j,k),v2c(m,l,i,j,k), - & v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k)) - enddo - enddo - enddo + read (itordp,*) ntermd_1(i,j,k,iblock), + & ntermd_2(i,j,k,iblock) + ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock) + ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock) + read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1, + & ntermd_1(i,j,k,iblock)) + read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1, + & ntermd_1(i,j,k,iblock)) + read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1, + & ntermd_1(i,j,k,iblock)) + read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1, + & ntermd_1(i,j,k,iblock)) + C Martix of D parameters for one dimesional foureir series + do l=1,ntermd_1(i,j,k,iblock) + v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock) + v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock) + v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock) + v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock) + c write(iout,*) "whcodze" , + c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock) + enddo + read (itordp,*) ((v2c(l,m,i,j,k,iblock), + & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock), + & v2s(m,l,i,j,k,iblock), + & m=1,l-1),l=1,ntermd_2(i,j,k,iblock)) + C Martix of D parameters for two dimesional fourier series + do l=1,ntermd_2(i,j,k,iblock) + do m=1,l-1 + v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock) + v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock) + v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock) + v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock) + enddo!m + enddo!l + enddo!k + enddo!j + enddo!i + enddo!iblock if (lprint) then write (iout,*) write (iout,*) 'Constants for double torsionals' - do i=1,ntortyp - do j=1,ntortyp - do k=1,ntortyp + do iblock=1,2 + do i=0,ntortyp-1 + do j=-ntortyp+1,ntortyp-1 + do k=-ntortyp+1,ntortyp-1 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k, - & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k) + & ' nsingle',ntermd_1(i,j,k,iblock), + & ' ndouble',ntermd_2(i,j,k,iblock) write (iout,*) write (iout,*) 'Single angles:' - do l=1,ntermd_1(i,j,k) + do l=1,ntermd_1(i,j,k,iblock) write (iout,'(i5,2f10.5,5x,2f10.5)') l, - & v1c(1,l,i,j,k),v1s(1,l,i,j,k), - & v1c(2,l,i,j,k),v1s(2,l,i,j,k) + & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock), + & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock) enddo write (iout,*) write (iout,*) 'Pairs of angles:' - write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k)) - do l=1,ntermd_2(i,j,k) + write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock)) + do l=1,ntermd_2(i,j,k,iblock) write (iout,'(i5,20f10.5)') - & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k)) + & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)) enddo write (iout,*) - write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k)) - do l=1,ntermd_2(i,j,k) + write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock)) + do l=1,ntermd_2(i,j,k,iblock) write (iout,'(i5,20f10.5)') - & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k)) + & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)) enddo write (iout,*) enddo enddo enddo + enddo endif #endif C Read of Side-chain backbone correlation parameters @@@ -589,6 -683,10 +683,10 @@@ CC C read (isccor,*) nsccortyp read (isccor,*) (isccortyp(i),i=1,ntyp) + do i=-ntyp,-1 + isccortyp(i)=-isccortyp(-i) + enddo + iscprol=isccortyp(20) c write (iout,*) 'ntortyp',ntortyp maxinter=3 cc maxinter is maximum interaction sites @@@ -597,12 -695,52 +695,52 @@@ do j=1,nsccortyp read (isccor,*) nterm_sccor(i,j),nlor_sccor(i,j) v0ijsccor=0.0d0 + v0ijsccor1=0.0d0 + v0ijsccor2=0.0d0 + v0ijsccor3=0.0d0 si=-1.0d0 - + nterm_sccor(-i,j)=nterm_sccor(i,j) + nterm_sccor(-i,-j)=nterm_sccor(i,j) + nterm_sccor(i,-j)=nterm_sccor(i,j) do k=1,nterm_sccor(i,j) read (isccor,*) kk,v1sccor(k,l,i,j) & ,v2sccor(k,l,i,j) + if (j.eq.iscprol) then + if (i.eq.isccortyp(10)) then + v1sccor(k,l,i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) + else + v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 + & +v2sccor(k,l,i,j)*dsqrt(0.75d0) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 + & +v1sccor(k,l,i,j)*dsqrt(0.75d0) + v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j) + v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j) + v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j) + endif + else + if (i.eq.isccortyp(10)) then + v1sccor(k,l,i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) + else + if (j.eq.isccortyp(10)) then + v1sccor(k,l,-i,j)=v1sccor(k,l,i,j) + v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j) + else + v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) + v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j) + v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j) + v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j) + endif + endif + endif v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j) + v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j) + v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j) + v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j) si=-si enddo do k=1,nlor_sccor(i,j) @@@ -611,7 -749,10 +749,10 @@@ v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &(1+vlor3sccor(k,i,j)**2) enddo - v0sccor(i,j)=v0ijsccor + v0sccor(l,i,j)=v0ijsccor + v0sccor(l,-i,j)=v0ijsccor1 + v0sccor(l,i,-j)=v0ijsccor2 + v0sccor(l,-i,-j)=v0ijsccor3 enddo enddo enddo @@@ -640,39 -781,96 +781,96 @@@ C 9/18/99 (AL) Read coefficients of th C interaction energy of the Gly, Ala, and Pro prototypes. C read (ifourier,*) nloctyp - do i=1,nloctyp + do i=0,nloctyp-1 read (ifourier,*) - read (ifourier,*) (b(ii,i),ii=1,13) + read (ifourier,*) (b(ii),ii=1,13) if (lprint) then write (iout,*) 'Type',i - write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13) + write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13) endif - B1(1,i) = b(3,i) - B1(2,i) = b(5,i) - B1tilde(1,i) = b(3,i) - B1tilde(2,i) =-b(5,i) - B2(1,i) = b(2,i) - B2(2,i) = b(4,i) - CC(1,1,i)= b(7,i) - CC(2,2,i)=-b(7,i) - CC(2,1,i)= b(9,i) - CC(1,2,i)= b(9,i) - Ctilde(1,1,i)=b(7,i) - Ctilde(1,2,i)=b(9,i) - Ctilde(2,1,i)=-b(9,i) - Ctilde(2,2,i)=b(7,i) - DD(1,1,i)= b(6,i) - DD(2,2,i)=-b(6,i) - DD(2,1,i)= b(8,i) - DD(1,2,i)= b(8,i) - Dtilde(1,1,i)=b(6,i) - Dtilde(1,2,i)=b(8,i) - Dtilde(2,1,i)=-b(8,i) - Dtilde(2,2,i)=b(6,i) - EE(1,1,i)= b(10,i)+b(11,i) - EE(2,2,i)=-b(10,i)+b(11,i) - EE(2,1,i)= b(12,i)-b(13,i) - EE(1,2,i)= b(12,i)+b(13,i) + B1(1,i) = b(3) + B1(2,i) = b(5) + B1(1,-i) = b(3) + B1(2,-i) = -b(5) + c b1(1,i)=0.0d0 + c b1(2,i)=0.0d0 + B1tilde(1,i) = b(3) + B1tilde(2,i) =-b(5) + B1tilde(1,-i) =-b(3) + B1tilde(2,-i) =b(5) + c b1tilde(1,i)=0.0d0 + c b1tilde(2,i)=0.0d0 + B2(1,i) = b(2) + B2(2,i) = b(4) + B2(1,-i) =b(2) + B2(2,-i) =-b(4) + + c b2(1,i)=0.0d0 + c b2(2,i)=0.0d0 + CC(1,1,i)= b(7) + CC(2,2,i)=-b(7) + CC(2,1,i)= b(9) + CC(1,2,i)= b(9) + CC(1,1,-i)= b(7) + CC(2,2,-i)=-b(7) + CC(2,1,-i)=-b(9) + CC(1,2,-i)=-b(9) + 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 + Ctilde(1,1,i)=b(7) + Ctilde(1,2,i)=b(9) + Ctilde(2,1,i)=-b(9) + Ctilde(2,2,i)=b(7) + Ctilde(1,1,-i)=b(7) + Ctilde(1,2,-i)=-b(9) + Ctilde(2,1,-i)=b(9) + Ctilde(2,2,-i)=b(7) + + 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 + DD(1,1,i)= b(6) + DD(2,2,i)=-b(6) + DD(2,1,i)= b(8) + DD(1,2,i)= b(8) + DD(1,1,-i)= b(6) + DD(2,2,-i)=-b(6) + DD(2,1,-i)=-b(8) + DD(1,2,-i)=-b(8) + 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 + Dtilde(1,1,i)=b(6) + Dtilde(1,2,i)=b(8) + Dtilde(2,1,i)=-b(8) + Dtilde(2,2,i)=b(6) + Dtilde(1,1,-i)=b(6) + Dtilde(1,2,-i)=-b(8) + Dtilde(2,1,-i)=b(8) + Dtilde(2,2,-i)=b(6) + + 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 + EE(1,1,i)= b(10)+b(11) + EE(2,2,i)=-b(10)+b(11) + EE(2,1,i)= b(12)-b(13) + EE(1,2,i)= b(12)+b(13) + EE(1,1,-i)= b(10)+b(11) + EE(2,2,-i)=-b(10)+b(11) + EE(2,1,-i)=-b(12)+b(13) + EE(1,2,-i)=-b(12)-b(13) + + 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 do i=1,nloctyp