From 2792b8ba89a535de351da0c8ef2f2fb1c491f42e Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Thu, 3 Aug 2017 14:17:37 +0200 Subject: [PATCH 1/1] split matrices working reading pdb --- source/unres/control.F90 | 88 ++++++++++++---------- source/unres/energy.f90 | 3 + source/unres/geometry.f90 | 3 +- source/unres/io.f90 | 3 +- source/unres/io_config.f90 | 176 ++++++++++++++++++++++++++++++++++++++++---- 5 files changed, 219 insertions(+), 54 deletions(-) diff --git a/source/unres/control.F90 b/source/unres/control.F90 index 87a0793..9b5f441 100644 --- a/source/unres/control.F90 +++ b/source/unres/control.F90 @@ -490,7 +490,7 @@ ind_scpint_old,nsumgrad,nlen,ngrad_start,ngrad_end,& ierror,k,ierr,iaux,ncheck_to,ncheck_from,ind_typ,& ichunk,int_index_old - + integer,dimension(5) :: nct_molec !el allocate(itask_cont_from(0:nfgtasks-1)) !(0:max_fg_procs-1) !el allocate(itask_cont_to(0:nfgtasks-1)) !(0:max_fg_procs-1) @@ -632,7 +632,11 @@ ispp=4 !?? wham ispp=2 #ifdef MPI ! Now partition the electrostatic-interaction array - npept=nct-nnt + if (itype(nres_molec(1),1).eq.ntyp1_molec(1)) then + npept=nres_molec(1)-nnt-1 + else + npept=nres_molec(1)-nnt + endif nele_int_tot=(npept-ispp)*(npept-ispp+1)/2 call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde) if (lprint) & @@ -644,7 +648,13 @@ iatel_e=0 ind_eleint=0 ind_eleint_old=0 - do i=nnt,nct-3 + if (itype(nres_molec(1),1).eq.ntyp_molec(1)) then + nct_molec(1)=nres_molec(1)-1 + else + nct_molec(1)=nres_molec(1) + endif + + do i=nnt,nct_molec(1)-3 ijunk=0 call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,& iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13) @@ -660,7 +670,7 @@ ind_eleint_vdw_old=0 iatel_s_vdw=0 iatel_e_vdw=0 - do i=nnt,nct-3 + do i=nnt,nct_molec(1)-3 ijunk=0 call int_partition(ind_eleint_vdw,my_ele_inds_vdw,& my_ele_inde_vdw,i,& @@ -673,16 +683,16 @@ 15 continue #else iatel_s=nnt - iatel_e=nct-5 ! ?? wham iatel_e=nct-3 + iatel_e=nct_molec(1)-5 ! ?? wham iatel_e=nct-3 do i=iatel_s,iatel_e ielstart(i)=i+4 ! ?? wham +2 - ielend(i)=nct-1 + ielend(i)=nct_molec(1)-1 enddo iatel_s_vdw=nnt - iatel_e_vdw=nct-3 + iatel_e_vdw=nct_molec(1)-3 do i=iatel_s_vdw,iatel_e_vdw ielstart_vdw(i)=i+2 - ielend_vdw(i)=nct-1 + ielend_vdw(i)=nct_molec(1)-1 enddo #endif if (lprint) then @@ -707,7 +717,7 @@ iatscp_e=0 ind_scpint=0 ind_scpint_old=0 - do i=nnt,nct-1 + do i=nnt,nct_molec(1)-1 if (i.lt.nnt+iscp) then !d write (iout,*) 'i.le.nnt+iscp' call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,& @@ -731,12 +741,12 @@ 14 continue #else iatscp_s=nnt - iatscp_e=nct-1 - do i=nnt,nct-1 + iatscp_e=nct_molec(1)-1 + do i=nnt,nct_molec(1)-1 if (i.lt.nnt+iscp) then nscp_gr(i)=1 iscpstart(i,1)=i+iscp - iscpend(i,1)=nct + iscpend(i,1)=nct_molec(1) elseif (i.gt.nct-iscp) then nscp_gr(i)=1 iscpstart(i,1)=nnt @@ -746,7 +756,7 @@ iscpstart(i,1)=nnt iscpend(i,1)=i-iscp iscpstart(i,2)=i+iscp - iscpend(i,2)=nct + iscpend(i,2)=nct_molec(1) endif enddo ! i #endif @@ -760,47 +770,48 @@ endif ! lprint ! Partition local interactions #ifdef MPI - call int_bounds(nres-2,loc_start,loc_end) + call int_bounds(nres_molec(1)-2,loc_start,loc_end) loc_start=loc_start+1 loc_end=loc_end+1 - call int_bounds(nres-2,ithet_start,ithet_end) + call int_bounds(nres_molec(1)-2,ithet_start,ithet_end) ithet_start=ithet_start+2 ithet_end=ithet_end+2 - call int_bounds(nct-nnt-2,iturn3_start,iturn3_end) + call int_bounds(nct_molec(1)-nnt-2,iturn3_start,iturn3_end) iturn3_start=iturn3_start+nnt iphi_start=iturn3_start+2 iturn3_end=iturn3_end+nnt iphi_end=iturn3_end+2 iturn3_start=iturn3_start-1 iturn3_end=iturn3_end-1 - call int_bounds(nres-3,itau_start,itau_end) + call int_bounds(nres_molec(1)-3,itau_start,itau_end) itau_start=itau_start+3 itau_end=itau_end+3 - call int_bounds(nres-3,iphi1_start,iphi1_end) + call int_bounds(nres_molec(1)-3,iphi1_start,iphi1_end) iphi1_start=iphi1_start+3 iphi1_end=iphi1_end+3 - call int_bounds(nct-nnt-3,iturn4_start,iturn4_end) + call int_bounds(nct_molec(1)-nnt-3,iturn4_start,iturn4_end) iturn4_start=iturn4_start+nnt iphid_start=iturn4_start+2 iturn4_end=iturn4_end+nnt iphid_end=iturn4_end+2 iturn4_start=iturn4_start-1 iturn4_end=iturn4_end-1 + print *,"TUTUTU",nres_molec(1),nres call int_bounds(nres-2,ibond_start,ibond_end) ibond_start=ibond_start+1 ibond_end=ibond_end+1 - call int_bounds(nct-nnt,ibondp_start,ibondp_end) + call int_bounds(nct_molec(1)-nnt,ibondp_start,ibondp_end) ibondp_start=ibondp_start+nnt ibondp_end=ibondp_end+nnt - call int_bounds1(nres-1,ivec_start,ivec_end) + call int_bounds1(nres_molec(1)-1,ivec_start,ivec_end) ! print *,"Processor",myrank,fg_rank,fg_rank1, ! & " ivec_start",ivec_start," ivec_end",ivec_end iset_start=loc_start+2 iset_end=loc_end+2 - call int_bounds(nres,ilip_start,ilip_end) + call int_bounds(nres_molec(1),ilip_start,ilip_end) ilip_start=ilip_start ilip_end=ilip_end - call int_bounds(nres-1,itube_start,itube_end) + call int_bounds(nres_molec(1)-1,itube_start,itube_end) itube_start=itube_start itube_end=itube_end if (ndih_constr.eq.0) then @@ -1309,17 +1320,17 @@ endif #else loc_start=2 - loc_end=nres-1 + loc_end=nres_molec(1)-1 ithet_start=3 - ithet_end=nres + ithet_end=nres_molec(1) iturn3_start=nnt - iturn3_end=nct-3 + iturn3_end=nct_molec(1)-3 iturn4_start=nnt - iturn4_end=nct-4 + iturn4_end=nct_molec(1)-4 iphi_start=nnt+3 - iphi_end=nct + iphi_end=nct_molec(1) iphi1_start=4 - iphi1_end=nres + iphi1_end=nres_molec(1) idihconstr_start=1 idihconstr_end=ndih_constr ithetaconstr_start=1 @@ -1327,21 +1338,21 @@ iphid_start=iphi_start iphid_end=iphi_end-1 itau_start=4 - itau_end=nres + itau_end=nres_molec(1) ibond_start=2 - ibond_end=nres-1 + ibond_end=nres_molec(1)-1 ibondp_start=nnt - ibondp_end=nct-1 + ibondp_end=nct_molec(1)-1 ivec_start=1 - ivec_end=nres-1 + ivec_end=nres_molec(1)-1 iset_start=3 - iset_end=nres+1 + iset_end=nres_molec(1)+1 iint_start=2 - iint_end=nres-1 + iint_end=nres_molec(1)-1 ilip_start=1 - ilip_end=nres + ilip_end=nres_molec(1) itube_start=1 - itube_end=nres + itube_end=nres_molec(1) #endif !el common /przechowalnia/ ! deallocate(iturn3_start_all) @@ -1734,7 +1745,8 @@ endif else if (molecule.eq.2) then do i=1,ntyp1_molec(molecule) - if (nam(1:1).eq.restyp(i,molecule)) then + print *,nam(1:1),restyp(i,molecule)(1:1) + if (nam(1:1).eq.restyp(i,molecule)(1:1)) then rescode=i return endif diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index f7020ff..c21d21d 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -5184,6 +5184,9 @@ usumsqder=usumsqder+ud(j)*uprod2 enddo estr=estr+uprod/usum + if (energy_dec) write (iout,*) & + "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,& + AKSC(1,iti),AKSC(1,iti)*diff*diff do j=1,3 gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres) enddo diff --git a/source/unres/geometry.f90 b/source/unres/geometry.f90 index 038dede..7b6c961 100644 --- a/source/unres/geometry.f90 +++ b/source/unres/geometry.f90 @@ -1,4 +1,4 @@ - module geometry + module geometry !----------------------------------------------------------------------------- use io_units use names @@ -3039,6 +3039,7 @@ integer :: i,j,ires,nscat real(kind=8),dimension(3,20) :: sccor real(kind=8) :: sccmj + print *,"I am in sccenter",ires,nscat do j=1,3 sccmj=0.0D0 do i=1,nscat diff --git a/source/unres/io.f90 b/source/unres/io.f90 index 8505d11..24549e8 100644 --- a/source/unres/io.f90 +++ b/source/unres/io.f90 @@ -742,6 +742,7 @@ allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres allocate(dc(3,0:2*maxres)) !(3,0:maxres2) allocate(itype(maxres,5)) !(maxres) + allocate(istype(maxres)) ! ! Zero out tables. ! @@ -996,7 +997,7 @@ do i=1,nres_molec(molec) istype(i)=sugarcode(sequence(i,molec)(1:1),i) - itype(i,1)=rescode(i,sequence(i,molec)(2:4),iscode,molec) + itype(i,1)=rescode(i,sequence(i,molec)(2:2),iscode,molec) enddo endif diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 277b6ba..cf23282 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -2314,7 +2314,7 @@ use control_data use compare_data use MPI_data - use control, only: rescode + use control, only: rescode,sugarcode ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.LOCAL' @@ -2327,7 +2327,7 @@ ! include 'COMMON.CONTROL' ! include 'COMMON.DISTFIT' ! include 'COMMON.SETUP' - integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,& + integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,& ! ishift_pdb logical :: lprn=.true.,fail real(kind=8),dimension(3) :: e1,e2,e3 @@ -2336,17 +2336,23 @@ character(len=5) :: atom character(len=80) :: card real(kind=8),dimension(3,20) :: sccor - integer :: kkk,lll,icha,kupa !rescode, + integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode, + integer :: isugar + character*1 :: sugar real(kind=8) :: cou + real(kind=8),dimension(3) :: ccc !el local varables integer,dimension(2,maxres/3) :: hfrag_alloc integer,dimension(4,maxres/3) :: bfrag_alloc real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm) - + real(kind=8),dimension(:,:), allocatable :: c_temporary + integer,dimension(:,:) , allocatable :: itype_temporary efree_temp=0.0d0 ibeg=1 ishift1=0 ishift=0 + molecule=1 + counter=0 ! write (2,*) "UNRES_PDB",unres_pdb ires=0 ires_old=0 @@ -2384,9 +2390,11 @@ else if (card(:3).eq.'TER') then ! End current chain ires_old=ires+2 + ishift=ishift+1 ishift1=ishift1+1 - itype(ires_old,1)=ntyp1 - itype(ires_old-1,1)=ntyp1 + itype(ires_old,molecule)=ntyp1_molec(molecule) + itype(ires_old-1,molecule)=ntyp1_molec(molecule) + nres_molec(molecule)=nres_molec(molecule)+2 ibeg=2 ! write (iout,*) "Chain ended",ires,ishift,ires_old if (unres_pdb) then @@ -2402,6 +2410,7 @@ ! Read free energy if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp ! Fish out the ATOM cards. +! write(iout,*) 'card',card(1:20) if (index(card(1:4),'ATOM').gt.0) then read (card(12:16),*) atom ! write (iout,*) "! ",atom," !",ires @@ -2436,6 +2445,7 @@ if (res.ne.'GLY' .and. res.ne. 'ACE') then ishift=ishift-1 itype(1,1)=ntyp1 + nres_molec(molecule)=nres_molec(molecule)+1 endif ires=ires-ishift+ishift1 ires_old=ires @@ -2448,17 +2458,27 @@ ishift1=ishift1-1 !!!!! ! write (iout,*) "New chain started",ires,ishift,ishift1,"!" ires=ires-ishift+ishift1 + print *,ires,ishift,ishift1 ires_old=ires ibeg=0 else ishift=ishift-(ires-ishift+ishift1-ires_old-1) ires=ires-ishift+ishift1 ires_old=ires - endif + endif + print *,'atom',ires,atom if (res.eq.'ACE' .or. res.eq.'NHE') then itype(ires,1)=10 else - itype(ires,1)=rescode(ires,res,0,1) + if (atom.eq.'CA '.or.atom.eq.'N ') then + molecule=1 + itype(ires,molecule)=rescode(ires,res,0,molecule) +! nres_molec(molecule)=nres_molec(molecule)+1 + else + molecule=2 + itype(ires,molecule)=rescode(ires,res(3:4),0,molecule) +! nres_molec(molecule)=nres_molec(molecule)+1 + endif endif else ires=ires-ishift+ishift1 @@ -2471,23 +2491,69 @@ if (atom.eq.'CA' .or. atom.eq.'CH3' .or. & res.eq.'NHE'.and.atom(:2).eq.'HN') then read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) +! print *,ires,ishift,ishift1 ! write (iout,*) "backbone ",atom #ifdef DEBUG write (iout,'(2i3,2x,a,3f8.3)') & ires,itype(ires,1),res,(c(j,ires),j=1,3) #endif iii=iii+1 + nres_molec(molecule)=nres_molec(molecule)+1 do j=1,3 sccor(j,iii)=c(j,ires) enddo + else if (.not.unres_pdb .and. (atom.eq."C1'" .or. & + atom.eq."C2'" .or. atom.eq."C3'" & + .or. atom.eq."C4'" .or. atom.eq."O4'")) then + read(card(31:54),'(3f8.3)') (ccc(j),j=1,3) +!c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3) + print *,ires,ishift,ishift1 + counter=counter+1 +! iii=iii+1 +! do j=1,3 +! sccor(j,iii)=c(j,ires) +! enddo + do j=1,3 + c(j,ires)=c(j,ires)+ccc(j)/5.0 + enddo + if (counter.eq.5) then +! iii=iii+1 + nres_molec(molecule)=nres_molec(molecule)+1 +! do j=1,3 +! sccor(j,iii)=c(j,ires) +! enddo + counter=0 + endif + print *, "ATOM",atom(1:3) + else if (atom(1:3).eq."C5'") then + read (card(19:19),'(a1)') sugar + isugar=sugarcode(sugar,ires) + if (ibeg.eq.1) then + istype(1)=isugar + else + istype(ires)=isugar + endif + if (unres_pdb) then + read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) + else + iii=iii+1 + read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) + endif ! write (*,*) card(23:27),ires,itype(ires,1) else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. & atom.ne.'N' .and. atom.ne.'C' .and. & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. & - atom.ne.'OXT' .and. atom(:2).ne.'3H') then + atom.ne.'OXT' .and. atom(:2).ne.'3H' & + .and. atom.ne.'P '.and. & + atom(1:1).ne.'H' .and. & + atom.ne.'OP1' .and. atom.ne.'OP2 ') then ! write (iout,*) "sidechain ",atom +! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3) + if ((molecule.ne.2).or.(atom(3:3).ne."'")) then +! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3) iii=iii+1 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) + endif endif endif enddo @@ -2496,6 +2562,11 @@ ! Calculate dummy residue coordinates inside the "chain" of a multichain ! system nres=ires + if (ires_old.ne.ires) nres_molec(molecule)=nres_molec(molecule)-2 +! print *,'I have', nres_molec(:) + + do k=1,5 + if (nres_molec(k).eq.0) cycle do i=2,nres-1 ! write (iout,*) i,itype(i,1) ! if (itype(i,1).eq.ntyp1) then @@ -2506,8 +2577,21 @@ ! dc(j,i)=c(j,i) ! enddo ! endif - if (itype(i,1).eq.ntyp1) then - if (itype(i+1,1).eq.ntyp1) then + if (itype(i,k).eq.ntyp1_molec(k)) then + if (itype(i+1,k).eq.ntyp1_molec(k)) then + if (itype(i+2,k).eq.0) then + print *,"masz sieczke" + do j=1,5 + if (itype(i+2,j).ne.0) then + itype(i+1,k)=0 + itype(i+1,j)=ntyp1_molec(j) + nres_molec(k)=nres_molec(k)-1 + nres_molec(j)=nres_molec(j)+1 + go to 3331 + endif + enddo + 3331 continue + endif ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false @@ -2556,6 +2640,7 @@ endif !itype.eq.ntyp1 enddo + enddo ! Calculate the CM of the last side chain. if (iii.gt.0) then if (unres_pdb) then @@ -2569,9 +2654,11 @@ ! nres=ires nsup=nres nstart_sup=1 + print *,"molecule",molecule if (itype(nres,1).ne.10) then nres=nres+1 - itype(nres,1)=ntyp1 + itype(nres,molecule)=ntyp1_molec(molecule) + nres_molec(molecule)=nres_molec(molecule)+1 if (unres_pdb) then ! 2/15/2013 by Adam: corrected insertion of the last dummy residue call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail) @@ -2591,6 +2678,8 @@ enddo endif endif + print *,'I have',nres, nres_molec(:) + !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb #ifdef WHAM_RUN if (nres.ne.nres0) then @@ -2638,6 +2727,7 @@ c(j,nres+1)=c(j,1) c(j,2*nres)=c(j,nres) enddo + if (itype(1,1).eq.ntyp1) then nsup=nsup-1 nstart_sup=2 @@ -2660,6 +2750,24 @@ enddo endif endif +! First lets assign correct dummy to correct type of chain +! 1) First residue + if (itype(1,1).eq.ntyp1) then + if (itype(2,1).eq.0) then + do j=2,5 + if (itype(2,j).ne.0) then + itype(1,1)=0 + itype(1,j)=ntyp1_molec(j) + nres_molec(1)=nres_molec(1)-1 + nres_molec(j)=nres_molec(j)+1 + go to 3231 + endif + enddo +3231 continue + endif + endif + print *,'I have',nres, nres_molec(:) + ! Copy the coordinates to reference coordinates ! do i=1,2*nres ! do j=1,3 @@ -2673,8 +2781,8 @@ write (iout,'(a,3(3x,a5),5x,3(3x,a5))') & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" do ires=1,nres - write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') & - restyp(itype(ires,1),1),ires,(c(j,ires),j=1,3),& + write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') & + (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),& (c(j,ires+nres),j=1,3) enddo endif @@ -2689,7 +2797,47 @@ (c(j,nres+ires),j=1,3) enddo endif +! NOW LETS ROCK! SORTING + allocate(c_temporary(3,2*nres)) + allocate(itype_temporary(nres,5)) + itype_temporary(:,:)=0 + seqalingbegin=1 + do k=1,5 + do i=1,nres + if (itype(i,k).ne.0) then + do j=1,3 + c_temporary(j,seqalingbegin)=c(j,i) + c_temporary(j,seqalingbegin+nres)=c(j,i+nres) + + enddo + itype_temporary(seqalingbegin,k)=itype(i,k) + seqalingbegin=seqalingbegin+1 + endif + enddo + enddo + do i=1,2*nres + do j=1,3 + c(j,i)=c_temporary(j,i) + enddo + enddo + do k=1,5 + do i=1,nres + itype(i,k)=itype_temporary(i,k) + enddo + enddo + if (lprn) then + write (iout,'(/a)') & + "Cartesian coordinates of the reference structure after sorting" + write (iout,'(a,3(3x,a5),5x,3(3x,a5))') & + "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + do ires=1,nres + write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') & + (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),& + (c(j,ires+nres),j=1,3) + enddo + endif + print *,seqalingbegin,nres if(.not.allocated(vbld)) then allocate(vbld(2*nres)) do i=1,2*nres -- 1.7.9.5