wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / source / unres / src_MD-M / csa.f
1       subroutine make_array
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'COMMON.IOUNITS'
5       include 'COMMON.CHAIN'
6       include 'COMMON.INTERACT'
7       include 'COMMON.CSA'
8
9 ccccccccccccccccccccccccc
10 c Level-2: group
11 ccccccccccccccccccccccccc
12
13       indg=0
14       do k=1,numch
15 ccccccccccccccccccccccccccccccccccccccccc
16 ! Groups the THETAs and the GAMMAs
17        do j=2,nres-1
18          indg=indg+1
19          if (j.lt.nres-1) then
20            ngroup(indg)=2
21          else
22            ngroup(indg)=1
23          endif
24          do i=1,ngroup(indg)
25           igroup(1,i,indg)=i
26           igroup(2,i,indg)=j
27           igroup(3,i,indg)=k
28          enddo
29        enddo
30 ccccccccccccccccccccccccccccccccccccccccc
31       enddo
32 ! Groups the ALPHAs and the BETAs
33       do k=1,numch
34        do j=2,nres-1
35         if(itype(j).ne.10) then
36          indg=indg+1
37          ngroup(indg)=2
38          do i=1,ngroup(indg)
39           igroup(1,i,indg)=i+2
40           igroup(2,i,indg)=j
41           igroup(3,i,indg)=k
42          enddo
43         endif
44        enddo
45       enddo
46
47       ntotgr=indg
48       write(iout,*) 
49       write(iout,*) "# of groups: ",ntotgr
50       do i=1,ntotgr
51        write(iout,41) i,ngroup(i),((igroup(k,j,i),k=1,3),j=1,ngroup(i))
52       enddo
53 !      close(iout)
54
55    40 format(i3,3x,3i3)
56    41 format(2i3,3x,6(3i3,2x))
57
58       return
59       end
60 ccccccccccccccccccccccccccccccccccccccccccccccccc
61 ccccccccccccccccccccccccccccccccccccccccccccccccc
62       subroutine make_ranvar(n,m,idum)
63       implicit real*8 (a-h,o-z)
64       include 'DIMENSIONS'
65       include 'COMMON.IOUNITS'
66       include 'COMMON.CHAIN'
67       include 'COMMON.VAR'
68       include 'COMMON.BANK'
69 c al      m=0
70       print *,'HOHOHOHO Make_RanVar!!!!!',n,m
71       itrial=0
72       do while(m.lt.n .and. itrial.le.10000)
73         itrial=itrial+1
74         jeden=1
75         call gen_rand_conf(jeden,*10)
76 !        call intout
77         m=m+1
78         do j=2,nres-1
79           dihang_in(1,j,1,m)=theta(j+1)
80           dihang_in(2,j,1,m)=phi(j+2)
81           dihang_in(3,j,1,m)=alph(j)
82           dihang_in(4,j,1,m)=omeg(j)
83         enddo  
84         dihang_in(2,nres-1,1,m)=0.0d0
85         goto 20
86   10    write (iout,*) 'Failed to generate conformation #',m+1,
87      &   ' itrial=',itrial
88   20    continue
89       enddo
90       print *,'Make_RanVar!!!!! m=',m,' itrial=',itrial
91
92       return
93       end
94 ccccccccccccccccccccccccccccccccccccccccccccccccc
95 ccccccccccccccccccccccccccccccccccccccccccccccccc
96       subroutine make_ranvar_reg(n,idum)
97       implicit real*8 (a-h,o-z)
98       include 'DIMENSIONS'
99       include 'COMMON.IOUNITS'
100       include 'COMMON.CHAIN'
101       include 'COMMON.VAR'
102       include 'COMMON.BANK'
103       include 'COMMON.GEO'
104       m=0
105       print *,'HOHOHOHO Make_RanVar!!!!!'
106       itrial=0
107       do while(m.lt.n .and. itrial.le.10000)
108         itrial=itrial+1
109         jeden=1
110         call gen_rand_conf(jeden,*10)
111 !        call intout
112         m=m+1
113         do j=2,nres-1
114           dihang_in(1,j,1,m)=theta(j+1)
115           dihang_in(2,j,1,m)=phi(j+2)
116           dihang_in(3,j,1,m)=alph(j)
117           dihang_in(4,j,1,m)=omeg(j)
118           if(m.le.n*0.1) then
119            dihang_in(1,j,1,m)=90.0*deg2rad
120            dihang_in(2,j,1,m)=50.0*deg2rad
121           endif
122         enddo  
123         dihang_in(2,nres-1,1,m)=0.0d0
124         goto 20
125   10    write (iout,*) 'Failed to generate conformation #',m+1,
126      &   ' itrial=',itrial
127   20    continue
128       enddo
129       print *,'Make_RanVar!!!!! m=',m,' itrial=',itrial
130
131       return
132       end
133 ccccccccccccccccccccccccccccccccccccccccccccccccc
134 ccccccccccccccccccccccccccccccccccccccccccccccccc
135       subroutine from_pdb(n,idum)
136 c This subroutine stores the UNRES int variables generated from 
137 c subroutine readpdb into the 1st conformation of in dihang_in.
138 c Subsequent n-1 conformations of dihang_in have identical values
139 c of theta and phi as the 1st conformation but random values for
140 c alph and omeg.
141 c The array cref (also generated from subroutine readpdb) is stored
142 c to crefjlee to be used for rmsd calculation in CSA, if necessary.
143       implicit real*8 (a-h,o-z)
144       include 'DIMENSIONS'
145       include 'COMMON.IOUNITS'
146       include 'COMMON.CHAIN'
147       include 'COMMON.VAR'
148       include 'COMMON.BANK'
149       include 'COMMON.GEO'
150
151       m=1
152       do j=2,nres-1
153         dihang_in(1,j,1,m)=theta(j+1)
154         dihang_in(2,j,1,m)=phi(j+2)
155         dihang_in(3,j,1,m)=alph(j)
156         dihang_in(4,j,1,m)=omeg(j)
157       enddo
158       dihang_in(2,nres-1,1,k)=0.0d0
159
160       do m=2,n
161        do k=2,nres-1
162         dihang_in(1,k,1,m)=dihang_in(1,k,1,1)
163         dihang_in(2,k,1,m)=dihang_in(2,k,1,1)
164         if(dabs(dihang_in(3,k,1,1)).gt.1.d-6) then
165          dihang_in(3,k,1,m)=90.d0*ran1(idum)+90.d0
166          dihang_in(3,k,1,m)=dihang_in(3,k,1,m)*deg2rad
167         endif
168         if(dabs(dihang_in(4,k,1,1)).gt.1.d-6) then
169          dihang_in(4,k,1,m)=360.d0*ran1(idum)-180.d0
170          dihang_in(4,k,1,m)=dihang_in(4,k,1,m)*deg2rad
171         endif
172        enddo
173       enddo
174
175 c Store cref to crefjlee (they are in COMMON.CHAIN).
176       do k=1,2*nres
177        do kk=1,3
178         kkk=1
179         crefjlee(kk,k)=cref(kk,k,kkk)
180        enddo
181       enddo
182
183       open(icsa_native_int,file=csa_native_int,status="old")
184       do m=1,n
185          write(icsa_native_int,*) m,e
186          write(icsa_native_int,200) 
187      &    (dihang_in(1,k,1,m)*rad2deg,k=2,nres-1)
188          write(icsa_native_int,200) 
189      &    (dihang_in(2,k,1,m)*rad2deg,k=2,nres-2)
190          write(icsa_native_int,200) 
191      &    (dihang_in(3,k,1,m)*rad2deg,k=2,nres-1)
192          write(icsa_native_int,200) 
193      &    (dihang_in(4,k,1,m)*rad2deg,k=2,nres-1)
194       enddo
195
196       do k=1,nres
197        write(icsa_native_int,200) (crefjlee(i,k),i=1,3)
198       enddo
199       close(icsa_native_int)
200
201   200 format (8f10.4)
202
203       return
204       end
205 ccccccccccccccccccccccccccccccccccccccccccccccccc
206 ccccccccccccccccccccccccccccccccccccccccccccccccc
207       subroutine from_int(n,mm,idum)
208       implicit real*8 (a-h,o-z)
209       include 'DIMENSIONS'
210       include 'COMMON.IOUNITS'
211       include 'COMMON.CHAIN'
212       include 'COMMON.VAR'
213       include 'COMMON.INTERACT'
214       include 'COMMON.BANK'
215       include 'COMMON.GEO'
216       include 'COMMON.CONTACTS'
217       integer ilen
218       external ilen
219       logical fail
220       double precision energia(0:n_ene)
221
222       open(icsa_native_int,file=csa_native_int,status="old")
223       read (icsa_native_int,*)
224       call read_angles(icsa_native_int,*10)
225       goto 11
226    10 write (iout,'(2a)') "CHUJ NASTAPIL - error in ",
227      &  csa_native_int(:ilen(csa_native_int))
228    11 continue
229       call intout
230       do j=2,nres-1
231         dihang_in(1,j,1,1)=theta(j+1)
232         dihang_in(2,j,1,1)=phi(j+2)
233         dihang_in(3,j,1,1)=alph(j)
234         dihang_in(4,j,1,1)=omeg(j)
235       enddo
236       dihang_in(2,nres-1,1,1)=0.0d0
237
238 c         read(icsa_native_int,*) ind,e
239 c         read(icsa_native_int,200) (dihang_in(1,k,1,1),k=2,nres-1)
240 c         read(icsa_native_int,200) (dihang_in(2,k,1,1),k=2,nres-2)
241 c         read(icsa_native_int,200) (dihang_in(3,k,1,1),k=2,nres-1)
242 c         read(icsa_native_int,200) (dihang_in(4,k,1,1),k=2,nres-1)
243 c         dihang_in(2,nres-1,1,1)=0.d0
244
245          maxsi=100
246          maxcount_fail=100
247
248          do m=mm+2,n
249 c          do k=2,nres-1
250 c           dihang_in(1,k,1,m)=dihang_in(1,k,1,1)
251 c           dihang_in(2,k,1,m)=dihang_in(2,k,1,1)
252 c           if(abs(dihang_in(3,k,1,1)).gt.1.d-3) then
253 c            dihang_in(3,k,1,m)=90.d0*ran1(idum)+90.d0
254 c           endif
255 c           if(abs(dihang_in(4,k,1,1)).gt.1.d-3) then
256 c            dihang_in(4,k,1,m)=360.d0*ran1(idum)-180.d0
257 c           endif
258 c          enddo
259 c           call intout
260            fail=.true.
261
262            icount_fail=0
263
264            DO WHILE (FAIL .AND. ICOUNT_FAIL .LE. MAXCOUNT_FAIL)
265
266            do i=nnt,nct
267              if (itype(i).ne.10) then
268 cd             print *,'i=',i,' itype=',itype(i),' theta=',theta(i+1)
269                fail=.true.
270                ii=0
271                do while (fail .and. ii .le. maxsi)
272                  call gen_side(itype(i),theta(i+1),alph(i),omeg(i),fail)
273                  ii = ii+1
274                enddo
275              endif
276            enddo
277            call chainbuild
278            call etotal(energia(0))
279            fail = (energia(0).ge.1.0d20)
280            icount_fail=icount_fail+1
281
282            ENDDO
283
284            if (icount_fail.gt.maxcount_fail) then
285              write (iout,*) 
286      &       'Failed to generate non-overlaping near-native conf.',
287      &       m
288            endif
289
290            do j=2,nres-1
291              dihang_in(1,j,1,m)=theta(j+1)
292              dihang_in(2,j,1,m)=phi(j+2)
293              dihang_in(3,j,1,m)=alph(j)
294              dihang_in(4,j,1,m)=omeg(j)
295            enddo
296            dihang_in(2,nres-1,1,m)=0.0d0
297          enddo
298
299 c      do m=1,n
300 c        write(icsa_native_int,*) m,e
301 c         write(icsa_native_int,200) (dihang_in(1,k,1,m),k=2,nres-1)
302 c         write(icsa_native_int,200) (dihang_in(2,k,1,m),k=2,nres-2)
303 c         write(icsa_native_int,200) (dihang_in(3,k,1,m),k=2,nres-1)
304 c         write(icsa_native_int,200) (dihang_in(4,k,1,m),k=2,nres-1)
305 c      enddo
306 c     close(icsa_native_int)
307
308 c      do m=mm+2,n
309 c       do i=1,4
310 c        do j=2,nres-1
311 c         dihang_in(i,j,1,m)=dihang_in(i,j,1,m)*deg2rad
312 c        enddo
313 c       enddo
314 c      enddo
315
316       call dihang_to_c(dihang_in(1,1,1,1))
317
318 c Store c to cref (they are in COMMON.CHAIN).
319       do k=1,2*nres
320        do kk=1,3
321         crefjlee(kk,k)=c(kk,k)
322        enddo
323       enddo
324
325       call contact(.true.,ncont_ref,icont_ref,co)
326
327 c      do k=1,nres
328 c       write(icsa_native_int,200) (crefjlee(i,k),i=1,3)
329 c      enddo
330       close(icsa_native_int)
331
332   200 format (8f10.4)
333
334       return
335       end
336 ccccccccccccccccccccccccccccccccccccccccccccccccc
337 ccccccccccccccccccccccccccccccccccccccccccccccccc
338       subroutine dihang_to_c(aarray)
339       implicit real*8 (a-h,o-z)
340       include 'DIMENSIONS'
341       include 'COMMON.CSA'
342       include 'COMMON.BANK'
343       include 'COMMON.CHAIN'
344       include 'COMMON.GEO'
345       include 'COMMON.VAR'
346
347       dimension aarray(mxang,maxres,mxch)
348
349 c     do i=4,nres
350 c      phi(i)=dihang_in(1,i-2,1,1)
351 c     enddo
352       do i=2,nres-1
353        theta(i+1)=aarray(1,i,1)
354        phi(i+2)=aarray(2,i,1)
355        alph(i)=aarray(3,i,1)
356        omeg(i)=aarray(4,i,1)
357       enddo
358
359       call chainbuild
360
361       return
362       end
363 ccccccccccccccccccccccccccccccccccccccccccccccccc
364 ccccccccccccccccccccccccccccccccccccccccccccccccc