From 2ec03757b68a8fb1a604a7f84f83086291fc75db Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Wed, 29 Jan 2014 21:48:46 +0100 Subject: [PATCH] DO zmiany czytanie pisanie i cutoff --- .gitignore | 1 + source/unres/src_MD-M/COMMON.CHAIN | 2 + source/unres/src_MD-M/COMMON.INTERACT | 6 +- source/unres/src_MD-M/energy_p_new_barrier.F | 400 +++++++++++++++++++++++--- source/unres/src_MD-M/parmread.F | 4 +- source/unres/src_MD-M/readpdb.F | 61 +++- source/unres/src_MD-M/readrtns_CSA.F | 6 + 7 files changed, 428 insertions(+), 52 deletions(-) diff --git a/.gitignore b/.gitignore index d5a5c53..77a69e0 100644 --- a/.gitignore +++ b/.gitignore @@ -27,3 +27,4 @@ bin/unres/MD/unres_ifort_MPICH_GAB_czyt.exe bin/unres/MD-M/unres_Tc_procor_newparm_em64-D-symetr.exe DIL/ compinfo +period_build diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN index f343887..bf5f688 100644 --- a/source/unres/src_MD-M/COMMON.CHAIN +++ b/source/unres/src_MD-M/COMMON.CHAIN @@ -14,3 +14,5 @@ & nsup,nstart_sup,nstart_seq, & chain_length,iprzes,tabperm(maxperm,maxsym),nperm common /from_zscore/ nz_start,nz_end,iz_sc + double precision boxxsize,boxysize,boxzsize + common /box/ boxxsize,boxysize,boxzsize diff --git a/source/unres/src_MD-M/COMMON.INTERACT b/source/unres/src_MD-M/COMMON.INTERACT index 83af3fb..37d3e88 100644 --- a/source/unres/src_MD-M/COMMON.INTERACT +++ b/source/unres/src_MD-M/COMMON.INTERACT @@ -24,8 +24,10 @@ C 12/1/95 Array EPS included in the COMMON block. & rpp(2,2),epp(2,2),elpp6(2,2),elpp3(2,2),eps_scp(ntyp,2), & rscp(ntyp,2) c 12/5/03 modified 09/18/03 Bond stretching parameters. - double precision vbldp0,vbldsc0,akp,aksc,abond0,distchainmax + double precision vbldp0,vbldpDUM, + & vbldsc0,akp,aksc,abond0,distchainmax integer nbondterm - common /stretch/ vbldp0,vbldsc0(maxbondterm,ntyp),akp, + common /stretch/ vbldp0,vbldpDUM, + & vbldsc0(maxbondterm,ntyp),akp, & aksc(maxbondterm,ntyp),abond0(maxbondterm,ntyp), & distchainmax,nbondterm(ntyp) diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index d4e2092..635a41e 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -1413,6 +1413,7 @@ C include 'COMMON.CALC' include 'COMMON.CONTROL' logical lprn + integer xshift,yshift,zshift evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -1420,6 +1421,11 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon lprn=.false. c if (icall.eq.0) lprn=.false. ind=0 +C the loop over all 27 posible neigbours (for xshift=0,yshift=0,zshift=0 +C we have the original box) + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 do i=iatsc_s,iatsc_e itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle @@ -1427,6 +1433,32 @@ c if (icall.eq.0) lprn=.false. xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) +C Return atom into box, boxxsize is size of box in x dimension + 134 continue + if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize + if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize +C Condition for being inside the proper box + if ((xi.gt.((xshift+0.5d0)*boxxsize)).or. + & (xi.lt.((xshift-0.5d0)*boxxsize))) then + go to 134 + endif + 135 continue + if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize + if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize +C Condition for being inside the proper box + if ((yi.gt.((yshift+0.5d0)*boxysize)).or. + & (yi.lt.((yshift-0.5d0)*boxysize))) then + go to 135 + endif + 136 continue + if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxxsize + if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxxsize +C Condition for being inside the proper box + if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. + & (zi.lt.((zshift-0.5d0)*boxzsize))) then + go to 136 + endif + dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -1467,12 +1499,41 @@ c chip12=0.0D0 c alf1=0.0D0 c alf2=0.0D0 c alf12=0.0D0 - xj=c(1,nres+j)-xi - yj=c(2,nres+j)-yi - zj=c(3,nres+j)-zi + xj=c(1,nres+j) + yj=c(2,nres+j) + zj=c(3,nres+j) +C Return atom J into box the original box + 137 continue + if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize + if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize +C Condition for being inside the proper box + if ((xj.gt.((0.5d0)*boxxsize)).or. + & (xj.lt.((-0.5d0)*boxxsize))) then + go to 137 + endif + 138 continue + if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize + if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize +C Condition for being inside the proper box + if ((yj.gt.((0.5d0)*boxysize)).or. + & (yj.lt.((-0.5d0)*boxysize))) then + go to 138 + endif + 139 continue + if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize + if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize +C Condition for being inside the proper box + if ((zj.gt.((0.5d0)*boxzsize)).or. + & (zj.lt.((-0.5d0)*boxzsize))) then + go to 139 + endif + dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) + xj=xj-xi + yj=yj-yi + zj=zj-zi c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi c write (iout,*) "j",j," dc_norm", c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j) @@ -1536,6 +1597,9 @@ C Calculate angular part of the gradient. enddo ! j enddo ! iint enddo ! i + enddo ! zshift + enddo ! yshift + enddo ! xshift c write (iout,*) "Number of loop steps in EGB:",ind cccc energy_dec=.false. return @@ -2788,6 +2852,7 @@ c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms C C Loop over i,i+2 and i,i+3 pairs of the peptide groups C +C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition do i=iturn3_start,iturn3_end if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle @@ -2800,6 +2865,31 @@ C xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi +C Return atom into box, boxxsize is size of box in x dimension + 184 continue + if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize + if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize +C Condition for being inside the proper box + if ((xmedi.gt.((0.5d0)*boxxsize)).or. + & (xmedi.lt.((-0.5d0)*boxxsize))) then + go to 184 + endif + 185 continue + if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize + if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize +C Condition for being inside the proper box + if ((ymedi.gt.((0.5d0)*boxysize)).or. + & (ymedi.lt.((-0.5d0)*boxysize))) then + go to 185 + endif + 186 continue + if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxxsize + if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxxsize +C Condition for being inside the proper box + if ((zmedi.gt.((0.5d0)*boxzsize)).or. + & (zmedi.lt.((-0.5d0)*boxzsize))) then + go to 186 + endif num_conti=0 call eelecij(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) @@ -2818,12 +2908,42 @@ C xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi +C Return atom into box, boxxsize is size of box in x dimension + 194 continue + if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize + if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize +C Condition for being inside the proper box + if ((xmedi.gt.((0.5d0)*boxxsize)).or. + & (xmedi.lt.((-0.5d0)*boxxsize))) then + go to 194 + endif + 195 continue + if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize + if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize +C Condition for being inside the proper box + if ((ymedi.gt.((0.5d0)*boxysize)).or. + & (ymedi.lt.((-0.5d0)*boxysize))) then + go to 195 + endif + 196 continue + if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxxsize + if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxxsize +C Condition for being inside the proper box + if ((zmedi.gt.((0.5d0)*boxzsize)).or. + & (zmedi.lt.((-0.5d0)*boxzsize))) then + go to 196 + endif + num_conti=num_cont_hb(i) call eelecij(i,i+3,ees,evdw1,eel_loc) if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) & call eturn4(i,eello_turn4) num_cont_hb(i)=num_conti enddo ! i +C Loop over all neighbouring boxes + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c @@ -2838,6 +2958,32 @@ c xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi +C Return atom into box, boxxsize is size of box in x dimension + 164 continue + if (xmedi.gt.((xshift+0.5d0)*boxxsize)) xmedi=xmedi-boxxsize + if (xmedi.lt.((xshift-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize +C Condition for being inside the proper box + if ((xmedi.gt.((xshift+0.5d0)*boxxsize)).or. + & (xmedi.lt.((xshift-0.5d0)*boxxsize))) then + go to 164 + endif + 165 continue + if (ymedi.gt.((yshift+0.5d0)*boxysize)) ymedi=ymedi-boxysize + if (ymedi.lt.((yshift-0.5d0)*boxysize)) ymedi=ymedi+boxysize +C Condition for being inside the proper box + if ((ymedi.gt.((yshift+0.5d0)*boxysize)).or. + & (ymedi.lt.((yshift-0.5d0)*boxysize))) then + go to 165 + endif + 166 continue + if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxxsize + if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxxsize +C Condition for being inside the proper box + if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or. + & (zmedi.lt.((zshift-0.5d0)*boxzsize))) then + go to 166 + endif + c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) do j=ielstart(i),ielend(i) @@ -2847,6 +2993,10 @@ c write (iout,*) i,j,itype(i),itype(j) enddo ! j num_cont_hb(i)=num_conti enddo ! i + enddo ! zshift + enddo ! yshift + enddo ! xshift + c write (iout,*) "Number of loop steps in EELEC:",ind cd do i=1,nres cd write (iout,'(i3,3f10.5,5x,3f10.5)') @@ -2911,9 +3061,41 @@ c ind=ind+1 dx_normj=dc_norm(1,j) dy_normj=dc_norm(2,j) dz_normj=dc_norm(3,j) - xj=c(1,j)+0.5D0*dxj-xmedi - yj=c(2,j)+0.5D0*dyj-ymedi - zj=c(3,j)+0.5D0*dzj-zmedi +C xj=c(1,j)+0.5D0*dxj-xmedi +C yj=c(2,j)+0.5D0*dyj-ymedi +C zj=c(3,j)+0.5D0*dzj-zmedi + xj=c(1,j)+0.5D0*dxj + yj=c(2,j)+0.5D0*dyj + zj=c(3,j)+0.5D0*dzj +C if ((i+3).lt.j) then !this condition keeps for turn3 and turn4 not subject to PBC + 174 continue + if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize + if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize +C Condition for being inside the proper box + if ((xj.gt.((0.5d0)*boxxsize)).or. + & (xj.lt.((-0.5d0)*boxxsize))) then + go to 174 + endif + 175 continue + if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize + if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize +C Condition for being inside the proper box + if ((yj.gt.((0.5d0)*boxysize)).or. + & (yj.lt.((-0.5d0)*boxysize))) then + go to 175 + endif + 176 continue + if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize + if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize +C Condition for being inside the proper box + if ((zj.gt.((0.5d0)*boxzsize)).or. + & (zj.lt.((-0.5d0)*boxzsize))) then + go to 176 + endif +C endif !endPBC condintion + xj=xj-xmedi + yj=yj-ymedi + zj=zj-zmedi rij=xj*xj+yj*yj+zj*zj rrmij=1.0D0/rij rij=dsqrt(rij) @@ -3278,7 +3460,9 @@ cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij -c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) + if (eel_loc_ij.ne.0) + & write (iout,'(a4,2i4,8f9.5)')'chuj', + & i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4) eel_loc=eel_loc+eel_loc_ij C Partial derivatives in virtual-bond dihedral angles gamma @@ -3667,8 +3851,9 @@ c write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3 call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) s3=0.5d0*(pizda(1,1)+pizda(2,2)) eello_turn4=eello_turn4-(s1+s2+s3) - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'eturn4',i,j,-(s1+s2+s3) +c write(iout,*)'chujOWO', auxvec(1),b1(1,iti2) + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,3f7.3)') + & 'eturn4',i,j,-(s1+s2+s3),s1,s2,s3 cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), cd & ' eello_turn4_num',8*eello_turn4_num C Derivatives in gamma(i) @@ -3838,13 +4023,40 @@ C r0_scp=4.5d0 cd print '(a)','Enter ESCP' cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 do i=iatscp_s,iatscp_e if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) zi=0.5D0*(c(3,i)+c(3,i+1)) - +C Return atom into box, boxxsize is size of box in x dimension + 134 continue + if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize + if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize +C Condition for being inside the proper box + if ((xi.gt.((xshift+0.5d0)*boxxsize)).or. + & (xi.lt.((xshift-0.5d0)*boxxsize))) then + go to 134 + endif + 135 continue + if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize + if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize +C Condition for being inside the proper box + if ((yi.gt.((yshift+0.5d0)*boxysize)).or. + & (yi.lt.((yshift-0.5d0)*boxysize))) then + go to 135 + endif + 136 continue + if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxxsize + if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxxsize +C Condition for being inside the proper box + if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. + & (zi.lt.((zshift-0.5d0)*boxzsize))) then + go to 136 + endif do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) @@ -3855,9 +4067,36 @@ c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi c zj=c(3,nres+j)-zi C Uncomment following three lines for Ca-p interactions - xj=c(1,j)-xi - yj=c(2,j)-yi - zj=c(3,j)-zi + xj=c(1,j) + yj=c(2,j) + zj=c(3,j) + 174 continue + if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize + if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize +C Condition for being inside the proper box + if ((xj.gt.((0.5d0)*boxxsize)).or. + & (xj.lt.((-0.5d0)*boxxsize))) then + go to 174 + endif + 175 continue + if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize + if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize +C Condition for being inside the proper box + if ((yj.gt.((0.5d0)*boxysize)).or. + & (yj.lt.((-0.5d0)*boxysize))) then + go to 175 + endif + 176 continue + if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize + if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize +C Condition for being inside the proper box + if ((zj.gt.((0.5d0)*boxzsize)).or. + & (zj.lt.((-0.5d0)*boxzsize))) then + go to 176 + endif + xj=xj-xi + yj=yj-yi + zj=zj-zi rij=xj*xj+yj*yj+zj*zj r0ij=r0_scp r0ijsq=r0ij*r0ij @@ -3909,6 +4148,9 @@ cgrad enddo enddo ! iint enddo ! i + enddo !zshift + enddo !yshift + enddo !xshift return end C----------------------------------------------------------------------------- @@ -3934,13 +4176,40 @@ C evdw2_14=0.0d0 cd print '(a)','Enter ESCP' cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e + do xshift=-1,1 + do yshift=-1,1 + do zshift=-1,1 do i=iatscp_s,iatscp_e if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) zi=0.5D0*(c(3,i)+c(3,i+1)) - +C Return atom into box, boxxsize is size of box in x dimension + 134 continue + if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize + if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize +C Condition for being inside the proper box + if ((xi.gt.((xshift+0.5d0)*boxxsize)).or. + & (xi.lt.((xshift-0.5d0)*boxxsize))) then + go to 134 + endif + 135 continue + if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize + if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize +C Condition for being inside the proper box + if ((yi.gt.((yshift+0.5d0)*boxysize)).or. + & (yi.lt.((yshift-0.5d0)*boxysize))) then + go to 135 + endif + 136 continue + if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxxsize + if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxxsize +C Condition for being inside the proper box + if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. + & (zi.lt.((zshift-0.5d0)*boxzsize))) then + go to 136 + endif do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) @@ -3951,9 +4220,36 @@ c xj=c(1,nres+j)-xi c yj=c(2,nres+j)-yi c zj=c(3,nres+j)-zi C Uncomment following three lines for Ca-p interactions - xj=c(1,j)-xi - yj=c(2,j)-yi - zj=c(3,j)-zi + xj=c(1,j) + yj=c(2,j) + zj=c(3,j) + 174 continue + if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize + if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize +C Condition for being inside the proper box + if ((xj.gt.((0.5d0)*boxxsize)).or. + & (xj.lt.((-0.5d0)*boxxsize))) then + go to 174 + endif + 175 continue + if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize + if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize +C Condition for being inside the proper box + if ((yj.gt.((0.5d0)*boxysize)).or. + & (yj.lt.((-0.5d0)*boxysize))) then + go to 175 + endif + 176 continue + if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize + if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize +C Condition for being inside the proper box + if ((zj.gt.((0.5d0)*boxzsize)).or. + & (zj.lt.((-0.5d0)*boxzsize))) then + go to 176 + endif + xj=xj-xi + yj=yj-yi + zj=zj-zi rrij=1.0D0/(xj*xj+yj*yj+zj*zj) fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) @@ -4010,6 +4306,9 @@ cgrad enddo enddo ! iint enddo ! i + enddo !zshift + enddo !yshift + enddo !xshift do i=1,nct do j=1,3 gvdwc_scp(j,i)=expon*gvdwc_scp(j,i) @@ -4218,24 +4517,31 @@ c estr=0.0d0 estr1=0.0d0 do i=ibondp_start,ibondp_end - if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then - estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) - do j=1,3 - gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) - & *dc(j,i-1)/vbld(i) - enddo - if (energy_dec) write(iout,*) - & "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax) - else + if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle +c estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) +c do j=1,3 +c gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) +c & *dc(j,i-1)/vbld(i) +c enddo +c if (energy_dec) write(iout,*) +c & "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax) +c else +C Checking if it involves dummy (NH3+ or COO-) group + if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then +C YES vbldpDUM is the equlibrium length of spring for Dummy atom + diff = vbld(i)-vbldpDUM + else +C NO vbldp0 is the equlibrium lenght of spring for peptide group diff = vbld(i)-vbldp0 - if (energy_dec) write (iout,*) + endif + if (energy_dec) write (iout,'(a7,i5,4f7.3)') & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff estr=estr+diff*diff do j=1,3 gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i) enddo c write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3) - endif +c endif enddo estr=0.5d0*AKP*estr+estr1 c @@ -4373,6 +4679,7 @@ C In following comments this theta will be referred to as t_c. bthetk=bthet(k,itype2,ichir21,ichir22) endif thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k) +c write(iout,*) 'chuj tu', y(k),z(k) enddo dthett=thet_pred_mean*ssd thet_pred_mean=thet_pred_mean*ss+a0thet(it) @@ -4409,8 +4716,8 @@ C Derivatives of the "mean" values in gamma1 and gamma2. & E_theta,E_tc) endif etheta=etheta+ethetai - if (energy_dec) write (iout,'(a6,i5,0pf7.3)') - & 'ebend',i,ethetai + if (energy_dec) write (iout,'(a6,i5,0pf7.3,f7.3,i5)') + & 'ebend',i,ethetai,theta(i),itype(i) if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett) @@ -4431,7 +4738,8 @@ C--------------------------------------------------------------------------- C Calculate the contributions to both Gaussian lobes. C 6/6/97 - Deform the Gaussians using the factor of 1/(1+time) C The "polynomial part" of the "standard deviation" of this part of -C the distribution. +C the distributioni. + write (iout,*) thetai,thet_pred_mean sig=polthet(3,it) do j=2,0,-1 sig=sig*thet_pred_mean+polthet(j,it) @@ -4461,6 +4769,7 @@ C Following variable is sigma(t_c)**(-2) delthe0=thetai-theta0i term1=-0.5D0*sigcsq*delthec*delthec term2=-0.5D0*sig0inv*delthe0*delthe0 +C write (iout,*)'term1',term1,term2,sigcsq,delthec,sig0inv,delthe0 C Following fuzzy logic is to avoid underflows in dexp and subsequent INFs and C NaNs in taking the logarithm. We extract the largest exponent which is added C to the energy (this being the log of the distribution) at the end of energy @@ -4488,6 +4797,7 @@ C Contribution of the bending energy from this theta is just the -log of C the sum of the contributions from the two lobes and the pre-exponential C factor. Simple enough, isn't it? ethetai=(-dlog(termexp)-termm+dlog(termpre)) +C write (iout,*) 'termexp',termexp,termm,termpre,i C NOW the derivatives!!! C 6/6/97 Take into account the deformation. E_theta=(delthec*sigcsq*term1 @@ -4554,7 +4864,7 @@ C logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 do i=ithet_start,ithet_end - if (itype(i-1).eq.ntyp1) cycle + if ((itype(i-1).eq.ntyp1)) cycle if (iabs(itype(i+1)).eq.20) iblock=2 if (iabs(itype(i+1)).ne.20) iblock=1 dethetai=0.0d0 @@ -5581,8 +5891,12 @@ C Set lprn=.true. for debugging c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end - if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 - & .or. itype(i).eq.ntyp1) cycle +C ANY TWO ARE DUMMY ATOMS in row CYCLE + if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or. + & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)) .or. + & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle +C For introducing the NH3+ and COO- group please check the etor_d for reference +C and guidance etors_ii=0.0D0 if (iabs(itype(i)).eq.20) then iblock=2 @@ -5683,8 +5997,11 @@ c lprn=.true. etors_d=0.0D0 c write(iout,*) "a tu??" do i=iphid_start,iphid_end - if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 - & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle +C ANY TWO ARE DUMMY ATOMS in row CYCLE + if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or. + & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)).or. + & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1)) .or. + & ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1))) cycle itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) @@ -5694,9 +6011,20 @@ c write(iout,*) "a tu??" gloci2=0.0D0 iblock=1 if (iabs(itype(i+1)).eq.20) iblock=2 +C ADASKO: WHEN PARAMETERS FOR THIS TYPE OF BLOCKING GROUP IS READY UNCOMMENT +C CHECK WEATHER THERE IS NECCESITY FOR iblock=3 for COO- +C if (itype(i+1).eq.ntyp1) iblock=3 +C The problem of NH3+ group can be resolved by adding new parameters please note if there +C IS or IS NOT need for this +C IF Yes uncomment below and add to parmread.F appropriate changes and to v1cij and so on +C is (itype(i-3).eq.ntyp1) ntblock=2 +C ntblock is N-terminal blocking group C Regular cosine and sine terms do j=1,ntermd_1(itori,itori1,itori2,iblock) +C Example of changes for NH3+ blocking group +C do j=1,ntermd_1(itori,itori1,itori2,iblock,ntblock) +C v1cij=v1c(1,j,itori,itori1,itori2,iblock,ntblock) v1cij=v1c(1,j,itori,itori1,itori2,iblock) v1sij=v1s(1,j,itori,itori1,itori2,iblock) v2cij=v1c(2,j,itori,itori1,itori2,iblock) diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 6ea0840..bd2165b 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -59,7 +59,7 @@ c Read the virtual-bond parameters, masses, and moments of inertia c and Stokes' radii of the peptide group and side chains c #ifdef CRYST_BOND - read (ibond,*) vbldp0,akp,mp,ip,pstok + read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok do i=1,ntyp nbondterm(i)=1 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i) @@ -71,7 +71,7 @@ c endif enddo #else - read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok + read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok do i=1,ntyp read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i), & j=1,nbondterm(i)),msc(i),isc(i),restok(i) diff --git a/source/unres/src_MD-M/readpdb.F b/source/unres/src_MD-M/readpdb.F index f2daadd..6b385fe 100644 --- a/source/unres/src_MD-M/readpdb.F +++ b/source/unres/src_MD-M/readpdb.F @@ -46,7 +46,8 @@ crc---------------------------------------- goto 10 else if (card(:3).eq.'TER') then C End current chain - ires_old=ires+1 + ires_old=ires+2 + itype(ires_old-1)=ntyp1 itype(ires_old)=ntyp1 ibeg=2 c write (iout,*) "Chain ended",ires,ishift,ires_old @@ -120,13 +121,49 @@ C system do i=2,nres-1 c write (iout,*) i,itype(i) if (itype(i).eq.ntyp1) then -c write (iout,*) "dummy",i,itype(i) - do j=1,3 - c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2 -c c(j,i)=(c(j,i-1)+c(j,i+1))/2 - dc(j,i)=c(j,i) - enddo - endif + if (itype(i+1).eq.ntyp1) then +C 16/01/2014 by Adasko: Adding to dummy atoms in the chain +C first is connected prevous chain (itype(i+1).eq.ntyp1)=true +C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false + if (unres_pdb) then +C 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 + 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 +C 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 + 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 C Calculate the CM of the last side chain. if (unres_pdb) then @@ -150,11 +187,11 @@ C 2/15/2013 by Adam: corrected insertion of the last dummy residue e2(3)=0.0d0 endif do j=1,3 - c(j,nres)=c(j,nres-1)-3.8d0*e2(j) + 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) + dcj=(c(j,nres-2)-c(j,nres-3))/2.0 c(j,nres)=c(j,nres-1)+dcj c(j,2*nres)=c(j,nres) enddo @@ -181,11 +218,11 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue e2(3)=0.0d0 endif do j=1,3 - c(j,1)=c(j,2)-3.8d0*e2(j) + c(j,1)=c(j,2)-1.9d0*e2(j) enddo else do j=1,3 - dcj=c(j,4)-c(j,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 diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 265b705..de4b92e 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -213,7 +213,13 @@ cfmc modecalc=10 i2ndstr=index(controlcard,'USE_SEC_PRED') gradout=index(controlcard,'GRADOUT').gt.0 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0 +C DISTCHAINMAX become obsolete for periodic boundry condition call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0) +C Reading the dimensions of box in x,y,z coordinates + call reada(controlcard,'BOXX',boxxsize,100.0d0) + call reada(controlcard,'BOXY',boxysize,100.0d0) + call reada(controlcard,'BOXZ',boxzsize,100.0d0) + if (me.eq.king .or. .not.out1file ) & write (iout,*) "DISTCHAINMAX",distchainmax -- 1.7.9.5