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