read2sigma in wham and cluster_wham
[unres.git] / source / wham / src-M / rmscalc.f
1       double precision function rmscalc(ishif,i,j,jcon,lprn)
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'DIMENSIONS.ZSCOPT'
5       include 'DIMENSIONS.COMPAR'
6       include 'COMMON.IOUNITS'
7       include 'COMMON.COMPAR'
8       include 'COMMON.CHAIN' 
9       include 'COMMON.INTERACT'
10       include 'COMMON.VAR'
11       include 'COMMON.CONTROL'
12       double precision przes(3),obrot(3,3)
13       double precision creff(3,maxres2),cc(3,maxres2)
14       logical iadded(maxres)
15       integer inumber(2,maxres)
16       common /ccc/ creff,cc,iadded,inumber
17       logical lprn
18       logical non_conv
19       integer ishif,i,j
20       if (lprn) then
21         write (iout,*) "i",i," j",j," jcont",jcon," ishif",ishif
22         write (iout,*) "npiece",npiece(j,i)
23         call flush(iout)
24       endif
25 c      write (iout,*) "symetr",symetr
26 c      call flush(iout)
27       nperm=1
28       do idup=1,symetr
29       nperm=nperm*idup
30       enddo
31 c      write (iout,*) "nperm",nperm
32 c      call flush(iout)
33       do kkk=1,nperm
34       idup=0
35       do l=1,nres
36         iadded(l)=.false.
37       enddo
38 c      write (iout,*) "kkk",kkk
39 c      call flush(iout)
40       do k=1,npiece(j,i)
41         if (i.eq.1) then
42           if (lprn) then
43             write (iout,*) "Level 1: j=",j,"k=",k," adding fragment",
44      &         ifrag(1,k,j),ifrag(2,k,j)
45             call flush(iout)
46           endif
47           call cprep(ifrag(1,k,j),ifrag(2,k,j),ishif,idup,kkk)
48 c          write (iout,*) "Exit cprep"
49 c          call flush(iout)
50 c          write (iout,*) "ii=",ii
51         else
52           kk = ipiece(k,j,i)
53 c          write (iout,*) "kk",kk," npiece",npiece(kk,1)
54           do l=1,npiece(kk,1)
55             if (lprn) then
56               write (iout,*) "Level",i,": j=",j,"k=",k," kk=",kk,
57      &          " l=",l," adding fragment",
58      &          ifrag(1,l,kk),ifrag(2,l,kk)
59               call flush(iout)
60             endif
61             call cprep(ifrag(1,l,kk),ifrag(2,l,kk),ishif,idup,kkk)
62 c            write (iout,*) "After cprep"
63 c            call flush(iout)
64           enddo 
65         endif
66       enddo
67       enddo
68       if (lprn) then
69         write (iout,*) "tuszukaj"
70         do kkk=1,nperm
71           do k=1,idup
72             write(iout,'(5i4,2(3f10.5,5x))') i,j,k,inumber(1,k),
73      &        inumber(2,k),(creff(l,k),l=1,3),(cc(l,k),l=1,3)
74           enddo
75         enddo
76         call flush(iout)
77       endif
78       rminrms=1.0d10
79       do kkk=1,nperm
80       call fitsq(rms,cc(1,1),creff(1,1),idup,przes,obrot,non_conv)
81       if (non_conv) then
82         print *,'Error: FITSQ non-convergent, jcon',jcon,i
83         rms = 1.0d10
84       else if (rms.lt.-1.0d-6) then 
85         print *,'Error: rms^2 = ',rms,jcon,i
86         rms = 1.0d10
87       else if (rms.ge.1.0d-6 .and. rms.lt.0) then
88         rms = 0.0d0
89       endif
90 c      write (iout,*) "rmsmin", rminrms, "rms", rms
91       if (rms.le.rminrms) rminrms=rms
92       enddo
93       rmscalc = dsqrt(rminrms)
94 c      write (iout, *) "analysys", rmscalc,anatemp
95       return
96       end
97 c-------------------------------------------------------------------------
98       subroutine cprep(if1,if2,ishif,idup,kwa)
99       implicit real*8 (a-h,o-z)
100       include 'DIMENSIONS'
101       include 'DIMENSIONS.ZSCOPT'
102       include 'DIMENSIONS.COMPAR'
103       include 'COMMON.CONTROL'
104       include 'COMMON.IOUNITS'
105       include 'COMMON.COMPAR'
106       include 'COMMON.CHAIN' 
107       include 'COMMON.INTERACT'
108       include 'COMMON.VAR'
109       double precision przes(3),obrot(3,3)
110       double precision creff(3,maxres2),cc(3,maxres2)
111       logical iadded(maxres)
112       integer inumber(2,maxres),iistrart,kwa,blar
113       common /ccc/ creff,cc,iadded,inumber
114 c      write (iout,*) "Calling cprep symetr",symetr," kwa",kwa
115       nperm=1
116       do blar=1,symetr
117       nperm=nperm*blar
118       enddo
119 c      write (iout,*) "nperm",nperm
120       kkk=kwa
121 c      ii=0
122       do l=if1,if2
123 c        write (iout,*) "l",l," iadded",iadded(l)
124 c        call flush(iout)
125         if (l+ishif.gt.1 .and. l+ishif.le.nres .and. .not.iadded(l)) 
126      &  then
127           idup=idup+1
128           iadded(l)=.true.
129           inumber(1,idup)=l
130           inumber(2,idup)=l+ishif
131           do m=1,3
132             creff(m,idup)=cref(m,l,kkk)
133             cc(m,idup)=c(m,l+ishif)
134           enddo
135         endif
136       enddo
137       return
138       end
139 c-------------------------------------------------------------------------
140       double precision function rmsnat(jcon)
141       implicit real*8 (a-h,o-z)
142       include 'DIMENSIONS'
143       include 'DIMENSIONS.ZSCOPT'
144       include 'DIMENSIONS.COMPAR'
145       include 'COMMON.IOUNITS'
146       include 'COMMON.COMPAR'
147       include 'COMMON.CHAIN' 
148       include 'COMMON.INTERACT'
149       include 'COMMON.VAR'
150       include 'COMMON.CONTROL'
151       double precision przes(3),obrot(3,3),cc(3,2*maxres),
152      &   ccref(3,2*maxres)
153       logical non_conv
154       integer ishif,i,j,resprzesun
155       rminrms=10.0d10
156       rmsminsing=10d10
157       nperm=1
158       do i=1,symetr
159        nperm=nperm*i
160       enddo
161       do kkk=1,nperm
162        nnsup=0
163        do i=1,nres
164         if (itype(i).ne.21) then
165           nnsup=nnsup+1
166           do j=1,3
167             cc(j,nnsup)=c(j,i)
168             ccref(j,nnsup)=cref(j,i,kkk)
169           enddo
170         endif
171        enddo
172        call fitsq(rms,cc(1,1),ccref(1,1),nnsup,przes,obrot,non_conv)
173        if (non_conv) then
174         print *,'Error: FITSQ non-convergent, jcon',jcon,i
175         rms=1.0d10
176        else if (rms.lt.-1.0d-6) then 
177         print *,'Error: rms^2 = ',rms,jcon,i
178         rms = 1.0d10
179        else if (rms.ge.1.0d-6 .and. rms.lt.0) then
180         rms=0.0d0
181        endif
182        if (rms.le.rminrms) rminrms=rms
183 c       write (iout,*) "kkk",kkk," rmsnat",rms , rminrms
184       enddo
185       rmsnat = dsqrt(rminrms)
186 C      write (iout,*)  "analysys",rmsnat, anatemp
187 C      liczenie rmsdla pojedynczego lancucha
188       return
189       end
190 c-----------------------------------------------------------------------------
191       double precision function gyrate(jcon)
192       implicit real*8 (a-h,o-z)
193       include 'DIMENSIONS'
194       include 'COMMON.INTERACT'
195       include 'COMMON.CHAIN'
196       double precision cen(3),rg
197
198       do j=1,3
199        cen(j)=0.0d0
200       enddo
201
202       do i=nnt,nct
203           do j=1,3
204             cen(j)=cen(j)+c(j,i)
205           enddo
206       enddo
207       do j=1,3
208             cen(j)=cen(j)/dble(nct-nnt+1)
209       enddo
210       rg = 0.0d0
211       do i = nnt, nct
212         do j=1,3
213          rg = rg + (c(j,i)-cen(j))**2
214         enddo
215       end do
216       gyrate = dsqrt(rg/dble(nct-nnt+1))
217       return
218       end