Adam's changes
[unres.git] / source / wham / src-M-SAXS.safe / oligomer.F
1       subroutine oligomer
2       implicit none
3       include "DIMENSIONS"
4       include "COMMON.CHAIN"
5       include "COMMON.INTERACT"
6       include "COMMON.IOUNITS"
7       integer nchain,i,ii,ipi,ipj,ipmin,j,jmin,k,ix,iy,iz,
8      &  ixmin,iymin,izmin
9       logical newchain 
10       integer ichain(2,20),iper(20),iaux
11       double precision dchain,dchainmin,cmchain(3,20)
12       nchain=1
13       newchain=.false.
14       ichain(1,nchain)=1
15       do i=2,nres
16         if (itype(i).eq.ntyp1) then
17           if (.not.newchain) then
18             ichain(2,nchain)=i
19             if (i.lt.nres) nchain=nchain+1
20             newchain=.true.
21           else 
22             newchain=.false.
23             ichain(1,nchain)=i
24           endif
25         endif   
26       enddo
27 #ifdef DEBUG
28       write (iout,*) "Chains"
29       do i=1,nchain
30         write (iout,*) i,ichain(1,i),ichain(2,i)
31       enddo
32 #endif
33       cmchain=0.0d0
34       do i=1,nchain
35         ii=0
36         do j=ichain(1,i),ichain(2,i)
37           if (itype(j).eq.ntyp1) cycle
38           ii=ii+1
39           do k=1,3
40             cmchain(k,i)=cmchain(k,i)+c(k,j)
41           enddo
42         enddo
43         do k=1,3
44           cmchain(k,i)=cmchain(k,i)/ii
45         enddo
46       enddo
47       do i=1,nchain
48         iper(i)=i
49       enddo
50       do i=1,nchain-1
51         dchainmin=1.0d10
52         do j=i+1,nchain
53           ipi=iper(i)
54           ipj=iper(j)
55           do ix=-1,1
56             do iy=-1,1
57               do iz=-1,1
58                 dchain=(cmchain(1,ipj)-cmchain(1,ipi)+ix*boxxsize)**2+
59      &                 (cmchain(2,ipj)-cmchain(2,ipi)+iy*boxysize)**2+
60      &                 (cmchain(3,ipj)-cmchain(3,ipi)+iz*boxzsize)**2
61 c                write (iout,*) "i",i," ipi",ipi," j",j," ipj",ipj," d",
62 c     &               dsqrt(dchain)," dmin",dsqrt(dchainmin)," jmin",jmin
63                 if (dchain.lt.dchainmin) then
64                   dchainmin=dchain
65                   ixmin=ix
66                   iymin=iy
67                   izmin=iz
68                   jmin=j
69                 endif
70               enddo
71             enddo
72           enddo
73         enddo
74         if (ixmin.eq.0 .and. iymin.eq.0 .and. izmin.eq.0) cycle
75         ipj=iper(jmin)
76         cmchain(1,ipj)=cmchain(1,ipj)+ixmin*boxxsize
77         cmchain(2,ipj)=cmchain(2,ipj)+iymin*boxysize
78         cmchain(3,ipj)=cmchain(3,ipj)+izmin*boxzsize
79         do k=ichain(1,ipj),ichain(2,ipj)
80           c(1,k)=c(1,k)+ixmin*boxxsize
81           c(2,k)=c(2,k)+iymin*boxysize
82           c(3,k)=c(3,k)+izmin*boxzsize
83           c(1,k+nres)=c(1,k+nres)+ixmin*boxxsize
84           c(2,k+nres)=c(2,k+nres)+iymin*boxysize
85           c(3,k+nres)=c(3,k+nres)+izmin*boxzsize
86         enddo
87 c        write (iout,*) "jmin",jmin," ipj",ipj,
88 c     &   " ixmin",ixmin," iymin",iymin," izmin",izmin
89         iaux=iper(i+1)
90         iper(i+1)=iper(jmin)
91         iper(jmin)=iaux
92       enddo
93       return
94       end