From: Adam Sieradzan Date: Mon, 17 Feb 2014 18:27:57 +0000 (+0100) Subject: The wrong one gradient compontent... X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=892c42fe1b098aac20f41e47faf97b1b2479b647 The wrong one gradient compontent... cleaning in progress by A.Liwo other gradient OK --- 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 8f89ae9..e56e104 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -2404,13 +2404,13 @@ c if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then if (i.gt. nnt+2 .and. i.lt.nct+2) then iti = itortyp(itype(i-2)) else - iti=ntortyp+1 + iti=ntortyp endif c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then iti1 = itortyp(itype(i-1)) else - iti1=ntortyp+1 + iti1=ntortyp endif cd write (iout,*) '*******i',i,' iti1',iti cd write (iout,*) 'b1',b1(:,iti) @@ -2451,10 +2451,10 @@ c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (itype(i-1).le.ntyp) then iti1 = itortyp(itype(i-1)) else - iti1=ntortyp+1 + iti1=ntortyp endif else - iti1=ntortyp+1 + iti1=ntortyp endif do k=1,2 mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1) @@ -2869,7 +2869,11 @@ 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 + & .or. itype(i+2).eq.ntyp1 + & .or. itype(i+3).eq.ntyp1 +c & .or. itype(i-1).eq.ntyp1 + & .or. itype(i+4).eq.ntyp1 + & ) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2912,7 +2916,9 @@ C Condition for being inside the proper box do i=iturn4_start,iturn4_end if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 & .or. itype(i+3).eq.ntyp1 - & .or. itype(i+4).eq.ntyp1) cycle + & .or. itype(i+4).eq.ntyp1 + & .or. itype(i+5).eq.ntyp1 + & ) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2962,7 +2968,9 @@ c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c do i=iatel_s,iatel_e - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 + & .or. itype(i+2).eq.ntyp1 + & ) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -3002,7 +3010,9 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) do j=ielstart(i),ielend(i) c write (iout,*) i,j,itype(i),itype(j) - if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle + if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1 + & .or.itype(j+2).eq.ntyp1 + &) cycle call eelecij(i,j,ees,evdw1,eel_loc) enddo ! j num_cont_hb(i)=num_conti @@ -3186,9 +3196,15 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo - ggg(1)=facvdw*xj - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj + if (sss.gt.0.0) then + ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj + else + ggg(1)=0.0 + ggg(2)=0.0 + ggg(3)=0.0 + endif c do k=1,3 c ghalf=0.5D0*ggg(k) c gvdwpp(k,i)=gvdwpp(k,i)+ghalf @@ -3212,7 +3228,7 @@ C MARYSIA facvdw=(ev1+evdwij)*sss facel=(el1+eesij) fac1=fac - fac=-3*rrmij*(facvdw+facvdw+facel)+sssgrad*rmij*evdwij + fac=-3*rrmij*(facvdw+facvdw+facel) erij(1)=xj*rmij erij(2)=yj*rmij erij(3)=zj*rmij @@ -3241,9 +3257,9 @@ cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo c 9/28/08 AL Gradient compotents will be summed only at the end - ggg(1)=facvdw*xj*sss - ggg(2)=facvdw*yj*sss - ggg(3)=facvdw*zj*sss + ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) @@ -3479,7 +3495,8 @@ cgrad endif C Contribution to the local-electrostatic energy coming from the i-j pair eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) & +a33*muij(4) -cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij +c write (iout,*) 'i',i,' j',j,itype(i),itype(j), +c & ' eel_loc_ij',eel_loc_ij if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij @@ -4652,9 +4669,8 @@ c time12=1.0d0 etheta=0.0D0 c write (*,'(a,i2)') 'EBEND ICG=',icg do i=ithet_start,ithet_end - print *,i,itype(i-1),itype(i),itype(i-2) - if (itype(i-1).eq.ntyp1) cycle - print *,'wchodze',itype(i-1) + if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 + & .or.itype(i).eq.ntyp1) cycle 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) @@ -4671,7 +4687,7 @@ C Zero the energy function and its derivative at 0 or pi. ichir22=isign(1,itype(i)) endif - if (i.gt.3 .and. itype(i-2).ne.ntyp1) then + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -4684,7 +4700,7 @@ C Zero the energy function and its derivative at 0 or pi. y(1)=0.0D0 y(2)=0.0D0 endif - if (i.lt.nres .and. itype(i).ne.ntyp1) then + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -4692,8 +4708,8 @@ C Zero the energy function and its derivative at 0 or pi. z(1)=cos(phii1) #else phii1=phi(i+1) - z(1)=dcos(phii1) #endif + z(1)=dcos(phii1) z(2)=dsin(phii1) else z(1)=0.0D0 @@ -4897,7 +4913,10 @@ C etheta=0.0D0 do i=ithet_start,ithet_end c print *,i,itype(i-1),itype(i),itype(i-2) - if ((itype(i-1).eq.ntyp1)) cycle + if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 + & .or.itype(i).eq.ntyp1) cycle +C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF + if (iabs(itype(i+1)).eq.20) iblock=2 if (iabs(itype(i+1)).ne.20) iblock=1 dethetai=0.0d0 @@ -4909,7 +4928,7 @@ c print *,i,itype(i-1),itype(i),itype(i-2) coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(i-2).ne.ntyp1) then + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -4930,7 +4949,7 @@ C propagation of chirality for glycine type sinph1(k)=0.0d0 enddo endif - if (i.lt.nres .and. itype(i).ne.ntyp1) then + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -5828,9 +5847,9 @@ c lprn=.true. do i=iphi_start,iphi_end etors_ii=0.0D0 if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 - & .or. itype(i).eq.ntyp1) cycle - itori=itortyp(itype(i-2)) - itori1=itortyp(itype(i-1)) + & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle + itori=itortyp(itype(i-2)) + itori1=itortyp(itype(i-1)) phii=phi(i) gloci=0.0D0 C Proline-Proline pair is a special case... @@ -5928,8 +5947,9 @@ C ANY TWO ARE DUMMY ATOMS in row CYCLE c if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or. c & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)) .or. c & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle - if ((itype(i-3).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. - & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) cycle + if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 + & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle +C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF C For introducing the NH3+ and COO- group please check the etor_d for reference C and guidance etors_ii=0.0D0 @@ -6033,10 +6053,14 @@ c lprn=.true. c write(iout,*) "a tu??" do i=iphid_start,iphid_end 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 +C if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or. +C & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)).or. +C & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1)) .or. +C & ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1))) cycle + if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or. + & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1).or. + & (itype(i+1).eq.ntyp1)) cycle +C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) @@ -6046,6 +6070,7 @@ C ANY TWO ARE DUMMY ATOMS in row CYCLE gloci2=0.0D0 iblock=1 if (iabs(itype(i+1)).eq.20) iblock=2 +C Iblock=2 Proline type 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 @@ -7174,7 +7199,7 @@ C--------------------------------------------------------------------------- if (j.lt.nres-1) then itj1 = itortyp(itype(j+1)) else - itj1=ntortyp+1 + itj1=ntortyp endif do iii=1,2 dipi(iii,1)=Ub2(iii,i) @@ -7264,14 +7289,14 @@ C parallel orientation of the two CA-CA-CA frames. if (i.gt.1) then iti=itortyp(itype(i)) else - iti=ntortyp+1 + iti=ntortyp endif itk1=itortyp(itype(k+1)) itj=itortyp(itype(j)) if (l.lt.nres-1) then itl1=itortyp(itype(l+1)) else - itl1=ntortyp+1 + itl1=ntortyp endif C A1 kernel(j+1) A2T cd do iii=1,2 @@ -7417,7 +7442,7 @@ C Antiparallel orientation of the two CA-CA-CA frames. if (i.gt.1) then iti=itortyp(itype(i)) else - iti=ntortyp+1 + iti=ntortyp endif itk1=itortyp(itype(k+1)) itl=itortyp(itype(l)) @@ -7425,7 +7450,7 @@ C Antiparallel orientation of the two CA-CA-CA frames. if (j.lt.nres-1) then itj1=itortyp(itype(j+1)) else - itj1=ntortyp+1 + itj1=ntortyp endif C A2 kernel(j-1)T A1T call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i), @@ -8579,14 +8604,14 @@ C energy moment and not to the cluster cumulant. if (j.lt.nres-1) then itj1=itortyp(itype(j+1)) else - itj1=ntortyp+1 + itj1=ntortyp endif itk=itortyp(itype(k)) itk1=itortyp(itype(k+1)) if (l.lt.nres-1) then itl1=itortyp(itype(l+1)) else - itl1=ntortyp+1 + itl1=ntortyp endif #ifdef MOMENT s1=dip(4,jj,i)*dip(4,kk,k) @@ -8698,19 +8723,19 @@ cd write (2,*) 'eello_graph4: wturn6',wturn6 if (j.lt.nres-1) then itj1=itortyp(itype(j+1)) else - itj1=ntortyp+1 + itj1=ntortyp endif itk=itortyp(itype(k)) if (k.lt.nres-1) then itk1=itortyp(itype(k+1)) else - itk1=ntortyp+1 + itk1=ntortyp endif itl=itortyp(itype(l)) if (l.lt.nres-1) then itl1=itortyp(itype(l+1)) else - itl1=ntortyp+1 + itl1=ntortyp endif cd write (2,*) 'eello6_graph4:','i',i,' j',j,' k',k,' l',l cd write (2,*) 'iti',iti,' itj',itj,' itj1',itj1,' itk',itk, diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F index 70e35ab..1a5e34a 100644 --- a/source/unres/src_MD-M/initialize_p.F +++ b/source/unres/src_MD-M/initialize_p.F @@ -245,6 +245,22 @@ C Initialize the bridge arrays ihpb(i)=0 jhpb(i)=0 enddo +C Initialize correlation arrays + do i=-maxtor,maxtor + do k=1,2 + b1(k,i)=0.0 + b2(k,i)=0.0 + b1tilde(k,i)=0.0 +c b2tilde(k,i)=0.0 + do j=1,2 + CC(j,k,i)=0.0 + Ctilde(j,k,i)=0.0 + DD(j,k,i)=0.0 + Dtilde(j,k,i)=0.0 + EE(j,k,i)=0.0 + enddo + enddo + enddo C C Initialize timing. C diff --git a/source/unres/src_MD-M/intcartderiv.F b/source/unres/src_MD-M/intcartderiv.F index d715c19..12d7dee 100644 --- a/source/unres/src_MD-M/intcartderiv.F +++ b/source/unres/src_MD-M/intcartderiv.F @@ -99,7 +99,8 @@ c conventional formulas around 0 and 180. #else do i=4,nres #endif -c if (itype(i-1).eq.21 .or. itype(i-2).eq.21 ) cycle +c if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 +c & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle c the conventional case sint=dsin(theta(i)) sint1=dsin(theta(i-1)) @@ -124,8 +125,8 @@ c Obtaining the gamma derivatives from sine derivative ctgt=cost/sint ctgt1=cost1/sint1 cosg_inv=1.0d0/cosg - if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then - dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1) +c if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then + dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1) & -(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2) dphi(j,1,i)=cosg_inv*dsinphi(j,1,i) dsinphi(j,2,i)= @@ -136,13 +137,13 @@ c Obtaining the gamma derivatives from sine derivative & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i) c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1) dphi(j,3,i)=cosg_inv*dsinphi(j,3,i) - endif +c endif c Bug fixed 3/24/05 (AL) enddo c Obtaining the gamma derivatives from cosine derivative else do j=1,3 - if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then +c if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3* & dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp* & dc_norm(j,i-3))/vbld(i-2) @@ -155,7 +156,7 @@ c Obtaining the gamma derivatives from cosine derivative & dcostheta(j,2,i)-fac0*(dc_norm(j,i-3)-scalp* & dc_norm(j,i-1))/vbld(i) dphi(j,3,i)=-1/sing*dcosphi(j,3,i) - endif +c endif enddo endif enddo diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 22236bc..bd2165b 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -591,7 +591,6 @@ C read (itorp,*,end=113,err=113) ntortyp,nterm_old if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp) - itortyp(ntyp1)=0 do i=1,ntortyp do j=1,ntortyp read (itorp,'(a)')