77fb71f03f207bfc93d0217c19102e6c2d28a5e9
[unres.git] / source / unres / src_MD / src / 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         crefjlee(kk,k)=cref(kk,k)
179        enddo
180       enddo
181
182       open(icsa_native_int,file=csa_native_int,status="old")
183       do m=1,n
184          write(icsa_native_int,*) m,e
185          write(icsa_native_int,200) 
186      &    (dihang_in(1,k,1,m)*rad2deg,k=2,nres-1)
187          write(icsa_native_int,200) 
188      &    (dihang_in(2,k,1,m)*rad2deg,k=2,nres-2)
189          write(icsa_native_int,200) 
190      &    (dihang_in(3,k,1,m)*rad2deg,k=2,nres-1)
191          write(icsa_native_int,200) 
192      &    (dihang_in(4,k,1,m)*rad2deg,k=2,nres-1)
193       enddo
194
195       do k=1,nres
196        write(icsa_native_int,200) (crefjlee(i,k),i=1,3)
197       enddo
198       close(icsa_native_int)
199
200   200 format (8f10.4)
201
202       return
203       end
204 ccccccccccccccccccccccccccccccccccccccccccccccccc
205 ccccccccccccccccccccccccccccccccccccccccccccccccc
206       subroutine from_int(n,mm,idum)
207       implicit real*8 (a-h,o-z)
208       include 'DIMENSIONS'
209       include 'COMMON.IOUNITS'
210       include 'COMMON.CHAIN'
211       include 'COMMON.VAR'
212       include 'COMMON.INTERACT'
213       include 'COMMON.BANK'
214       include 'COMMON.GEO'
215       include 'COMMON.CONTACTS'
216       integer ilen
217       external ilen
218       logical fail
219       double precision energia(0:n_ene)
220
221       open(icsa_native_int,file=csa_native_int,status="old")
222       read (icsa_native_int,*)
223       call read_angles(icsa_native_int,*10)
224       goto 11
225    10 write (iout,'(2a)') "CHUJ NASTAPIL - error in ",
226      &  csa_native_int(:ilen(csa_native_int))
227    11 continue
228       call intout
229       do j=2,nres-1
230         dihang_in(1,j,1,1)=theta(j+1)
231         dihang_in(2,j,1,1)=phi(j+2)
232         dihang_in(3,j,1,1)=alph(j)
233         dihang_in(4,j,1,1)=omeg(j)
234       enddo
235       dihang_in(2,nres-1,1,1)=0.0d0
236
237 c         read(icsa_native_int,*) ind,e
238 c         read(icsa_native_int,200) (dihang_in(1,k,1,1),k=2,nres-1)
239 c         read(icsa_native_int,200) (dihang_in(2,k,1,1),k=2,nres-2)
240 c         read(icsa_native_int,200) (dihang_in(3,k,1,1),k=2,nres-1)
241 c         read(icsa_native_int,200) (dihang_in(4,k,1,1),k=2,nres-1)
242 c         dihang_in(2,nres-1,1,1)=0.d0
243
244          maxsi=100
245          maxcount_fail=100
246
247          do m=mm+2,n
248 c          do k=2,nres-1
249 c           dihang_in(1,k,1,m)=dihang_in(1,k,1,1)
250 c           dihang_in(2,k,1,m)=dihang_in(2,k,1,1)
251 c           if(abs(dihang_in(3,k,1,1)).gt.1.d-3) then
252 c            dihang_in(3,k,1,m)=90.d0*ran1(idum)+90.d0
253 c           endif
254 c           if(abs(dihang_in(4,k,1,1)).gt.1.d-3) then
255 c            dihang_in(4,k,1,m)=360.d0*ran1(idum)-180.d0
256 c           endif
257 c          enddo
258 c           call intout
259            fail=.true.
260
261            icount_fail=0
262
263            DO WHILE (FAIL .AND. ICOUNT_FAIL .LE. MAXCOUNT_FAIL)
264
265            do i=nnt,nct
266              if (itype(i).ne.10) then
267 cd             print *,'i=',i,' itype=',itype(i),' theta=',theta(i+1)
268                fail=.true.
269                ii=0
270                do while (fail .and. ii .le. maxsi)
271                  call gen_side(itype(i),theta(i+1),alph(i),omeg(i),fail)
272                  ii = ii+1
273                enddo
274              endif
275            enddo
276            call chainbuild
277            call etotal(energia(0))
278            fail = (energia(0).ge.1.0d20)
279            icount_fail=icount_fail+1
280
281            ENDDO
282
283            if (icount_fail.gt.maxcount_fail) then
284              write (iout,*) 
285      &       'Failed to generate non-overlaping near-native conf.',
286      &       m
287            endif
288
289            do j=2,nres-1
290              dihang_in(1,j,1,m)=theta(j+1)
291              dihang_in(2,j,1,m)=phi(j+2)
292              dihang_in(3,j,1,m)=alph(j)
293              dihang_in(4,j,1,m)=omeg(j)
294            enddo
295            dihang_in(2,nres-1,1,m)=0.0d0
296          enddo
297
298 c      do m=1,n
299 c        write(icsa_native_int,*) m,e
300 c         write(icsa_native_int,200) (dihang_in(1,k,1,m),k=2,nres-1)
301 c         write(icsa_native_int,200) (dihang_in(2,k,1,m),k=2,nres-2)
302 c         write(icsa_native_int,200) (dihang_in(3,k,1,m),k=2,nres-1)
303 c         write(icsa_native_int,200) (dihang_in(4,k,1,m),k=2,nres-1)
304 c      enddo
305 c     close(icsa_native_int)
306
307 c      do m=mm+2,n
308 c       do i=1,4
309 c        do j=2,nres-1
310 c         dihang_in(i,j,1,m)=dihang_in(i,j,1,m)*deg2rad
311 c        enddo
312 c       enddo
313 c      enddo
314
315       call dihang_to_c(dihang_in(1,1,1,1))
316
317 c Store c to cref (they are in COMMON.CHAIN).
318       do k=1,2*nres
319        do kk=1,3
320         crefjlee(kk,k)=c(kk,k)
321        enddo
322       enddo
323
324       call contact(.true.,ncont_ref,icont_ref,co)
325
326 c      do k=1,nres
327 c       write(icsa_native_int,200) (crefjlee(i,k),i=1,3)
328 c      enddo
329       close(icsa_native_int)
330
331   200 format (8f10.4)
332
333       return
334       end
335 ccccccccccccccccccccccccccccccccccccccccccccccccc
336 ccccccccccccccccccccccccccccccccccccccccccccccccc
337       subroutine dihang_to_c(aarray)
338       implicit real*8 (a-h,o-z)
339       include 'DIMENSIONS'
340       include 'COMMON.CSA'
341       include 'COMMON.BANK'
342       include 'COMMON.CHAIN'
343       include 'COMMON.GEO'
344       include 'COMMON.VAR'
345
346       dimension aarray(mxang,maxres,mxch)
347
348 c     do i=4,nres
349 c      phi(i)=dihang_in(1,i-2,1,1)
350 c     enddo
351       do i=2,nres-1
352        theta(i+1)=aarray(1,i,1)
353        phi(i+2)=aarray(2,i,1)
354        alph(i)=aarray(3,i,1)
355        omeg(i)=aarray(4,i,1)
356       enddo
357
358       call chainbuild
359
360       return
361       end
362 ccccccccccccccccccccccccccccccccccccccccccccccccc
363 ccccccccccccccccccccccccccccccccccccccccccccccccc