unres src_MD-M from branch Multichain
[unres.git] / source / cluster / unres / src / readrtns.F
1       subroutine read_control
2 C
3 C Read molecular data
4 C
5       include 'DIMENSIONS'
6       include 'sizesclu.dat'
7       include 'COMMON.IOUNITS'
8       include 'COMMON.TIME1'
9       include 'COMMON.THREAD'
10       include 'COMMON.SBRIDGE'
11       include 'COMMON.CONTROL'
12       include 'COMMON.CLUSTER'
13       include 'COMMON.CHAIN'
14       include 'COMMON.HEADER'
15       character*320 controlcard,ucase
16 #ifdef MPL
17       include 'COMMON.INFO'
18 #endif
19
20       read (INP,'(a80)') titel
21       call card_concat(controlcard)
22
23 C Set up the time limit (caution! The time must be input in minutes!)
24       call reada(controlcard,'TIMLIM',timlim,60.0D0) 
25       write (iout,'(a,f10.1)') 'Time limit (min):',timlim
26       timlim=60.0D0*timlim
27 #ifdef MPL
28       timlim=timlim/nprocs
29 #endif
30         
31       modecalc=0
32       call readi(controlcard,'NRES',nres,0)
33       call reada(controlcard,'TEMP',betaT,298.0d0)
34       write (iout,*) "temperature",betaT
35       betaT=1.0d0/(betaT*1.987D-3)
36       write (iout,*) "betaT",betaT
37       call readi(controlcard,'PDBOUT',outpdb,0)
38       call readi(controlcard,'MOL2OUT',outmol2,0)
39       refstr=(index(controlcard,'REFSTR').gt.0)
40       pdbref=(index(controlcard,'PDBREF').gt.0)
41       iscode=index(controlcard,'ONE_LETTER')
42       tree=(index(controlcard,'MAKE_TREE').gt.0)
43       min_var=(index(controlcard,'MINVAR').gt.0)
44       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
45       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
46       call readi(controlcard,'NCUT',ncut,1)
47       call reada(controlcard,'ECUT',ecut,10.0d0)
48       lgrp=(index(controlcard,'LGRP').gt.0)
49       caonly=(index(controlcard,'CA_ONLY').gt.0)
50       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
51       print *,"print_dist",print_dist
52       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
53       call readi(controlcard,'IOPT',iopt,2) 
54       print *,"iopt",iopt
55       lside = index(controlcard,"SIDE").gt.0
56       efree = index(controlcard,"EFREE").gt.0
57       print *,controlcard
58       call readi(controlcard,'ISTART',is,1)
59       print *,"is",is
60       call readi(controlcard,'IEND',ie,10000000)
61       print *,"ie",ie
62       call readi(controlcard,'ISAMPL',isampl,1)
63       print *,"isampl",isampl
64       call reada(controlcard,'TS',ts,1.0d10)
65       call reada(controlcard,'TE',te,1.0d10)
66       print *,"is",is," ie",ie
67       if (min_var) iopt=1
68       return
69       end
70 c--------------------------------------------------------------------------
71       subroutine molread
72 C
73 C Read molecular data.
74 C
75       include 'DIMENSIONS'
76       include 'COMMON.IOUNITS'
77       include 'COMMON.GEO'
78       include 'COMMON.VAR'
79       include 'COMMON.INTERACT'
80       include 'COMMON.LOCAL'
81       include 'COMMON.NAMES'
82       include 'COMMON.CHAIN'
83       include 'COMMON.FFIELD'
84       include 'COMMON.SBRIDGE'
85       include 'COMMON.HEADER'
86       include 'COMMON.CONTROL'
87       include 'COMMON.DBASE'
88       include 'COMMON.THREAD'
89       include 'COMMON.CONTACTS'
90       include 'COMMON.TIME1'
91 #ifdef MPL
92       include 'COMMON.INFO'
93 #endif
94       character*4 sequence(maxres)
95       integer rescode
96       double precision x(maxvar)
97       character*64 pdbfile
98       dimension itype_pdb(maxres)
99       logical seq_comp
100 C
101 C Body
102 C
103       if (indpdb.gt.0 .or. pdbref) then
104         read(inp,'(a)') pdbfile
105         write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
106         open(ipdbin,file=pdbfile,status='old',err=33)
107         goto 34 
108   33    write (iout,'(a)') 'Error opening PDB file.'
109         stop
110   34    continue
111         print *,'Begin reading pdb data'
112         call readpdb
113         print *,'Finished reading pdb data'
114         write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
115         do i=1,nres
116           itype_pdb(i)=itype(i)
117         enddo
118         close (ipdbin)
119       endif
120       if (indpdb.eq.0) then
121 C Read sequence if not taken from the pdb file.
122         if (iscode.gt.0) then
123           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
124         else
125           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
126         endif
127 C Convert sequence to numeric code
128         do i=1,nres
129           itype(i)=rescode(i,sequence(i),iscode)
130         enddo
131       endif
132       print *,nres
133       print '(20i4)',(itype(i),i=1,nres)
134       do i=1,nres
135         if (itype(i).ne.20) then
136           itel(i)=1
137         else
138           itel(i)=2
139         endif  
140       enddo
141       print *,'Call Read_Bridge.'
142       call read_bridge
143       nnt=1
144       nct=nres
145       print *,'NNT=',NNT,' NCT=',NCT
146       if (itype(1).eq.21) nnt=2
147       if (itype(nres).eq.21) nct=nct-1
148       if (pdbref) then
149         write (iout,'(a,i3)') 'nsup=',nsup
150         nstart_seq=nnt
151         if (nsup.le.(nct-nnt+1)) then
152           do i=0,nct-nnt+1-nsup
153             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
154               nstart_seq=nnt+i
155               goto 111
156             endif
157           enddo
158           write (iout,'(a)') 
159      &            'Error - sequences to be superposed do not match.'
160           stop
161         else
162           do i=0,nsup-(nct-nnt+1)
163             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) 
164      &      then
165               nstart_sup=nstart_sup+i
166               nsup=nct-nnt+1
167               goto 111
168             endif
169           enddo 
170           write (iout,'(a)') 
171      &            'Error - sequences to be superposed do not match.'
172         endif
173   111   continue
174         write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
175      &                 ' nstart_seq=',nstart_seq
176       endif
177       if (refstr) then
178         if (.not.pdbref) then
179           call read_angles(inp,*38)
180           goto 39
181    38     write (iout,'(a)') 'Error reading reference structure.'
182 #ifdef MPL
183           call mp_stopall(Error_Msg)
184 #else
185           stop 'Error reading reference structure'
186 #endif
187    39     call chainbuild     
188           nstart_sup=nnt
189           nstart_seq=nnt
190           nsup=nct-nnt+1
191           do i=1,2*nres
192             do j=1,3
193               cref(j,i)=c(j,i)
194             enddo
195           enddo
196         endif
197         call contact(.true.,ncont_ref,icont_ref)
198       endif
199       return
200       end
201 c-----------------------------------------------------------------------------
202       logical function seq_comp(itypea,itypeb,length)
203       implicit none
204       integer length,itypea(length),itypeb(length)
205       integer i
206       do i=1,length
207         if (itypea(i).ne.itypeb(i)) then
208           seq_comp=.false.
209           return
210         endif
211       enddo
212       seq_comp=.true.
213       return
214       end
215 c-----------------------------------------------------------------------------
216       subroutine read_bridge
217 C Read information about disulfide bridges.
218       include 'DIMENSIONS'
219       include 'COMMON.IOUNITS'
220       include 'COMMON.GEO'
221       include 'COMMON.VAR'
222       include 'COMMON.INTERACT'
223       include 'COMMON.LOCAL'
224       include 'COMMON.NAMES'
225       include 'COMMON.CHAIN'
226       include 'COMMON.FFIELD'
227       include 'COMMON.SBRIDGE'
228       include 'COMMON.HEADER'
229       include 'COMMON.CONTROL'
230       include 'COMMON.DBASE'
231       include 'COMMON.THREAD'
232       include 'COMMON.TIME1'
233 #ifdef MPL
234       include 'COMMON.INFO'
235 #endif
236 C Read bridging residues.
237       read (inp,*) ns,(iss(i),i=1,ns)
238       print *,'ns=',ns
239 C Check whether the specified bridging residues are cystines.
240       do i=1,ns
241         if (itype(iss(i)).ne.1) then
242           write (iout,'(2a,i3,a)') 
243      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
244      &   ' can form a disulfide bridge?!!!'
245           write (*,'(2a,i3,a)') 
246      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
247      &   ' can form a disulfide bridge?!!!'
248 #ifdef MPL
249          call mp_stopall(error_msg)
250 #else
251          stop
252 #endif
253         endif
254       enddo
255 C Read preformed bridges.
256       if (ns.gt.0) then
257       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
258       if (nss.gt.0) then
259         nhpb=nss
260 C Check if the residues involved in bridges are in the specified list of
261 C bridging residues.
262         do i=1,nss
263           do j=1,i-1
264             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
265      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
266               write (iout,'(a,i3,a)') 'Disulfide pair',i,
267      &      ' contains residues present in other pairs.'
268               write (*,'(a,i3,a)') 'Disulfide pair',i,
269      &      ' contains residues present in other pairs.'
270 #ifdef MPL
271               call mp_stopall(error_msg)
272 #else
273               stop 
274 #endif
275             endif
276           enddo
277           do j=1,ns
278             if (ihpb(i).eq.iss(j)) goto 10
279           enddo
280           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
281    10     continue
282           do j=1,ns
283             if (jhpb(i).eq.iss(j)) goto 20
284           enddo
285           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
286    20     continue
287           dhpb(i)=dbr
288           forcon(i)=fbr
289         enddo
290         do i=1,nss
291           ihpb(i)=ihpb(i)+nres
292           jhpb(i)=jhpb(i)+nres
293         enddo
294       endif
295       endif
296       return
297       end
298 c----------------------------------------------------------------------------
299       subroutine read_angles(kanal,*)
300       include 'DIMENSIONS'
301       include 'COMMON.GEO'
302       include 'COMMON.VAR'
303       include 'COMMON.CHAIN'
304       include 'COMMON.IOUNITS'
305       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
306       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
307       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
308       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
309       do i=1,nres
310         theta(i)=deg2rad*theta(i)
311         phi(i)=deg2rad*phi(i)
312         alph(i)=deg2rad*alph(i)
313         omeg(i)=deg2rad*omeg(i)
314       enddo
315       return
316    10 return1
317       end
318 c----------------------------------------------------------------------------
319       subroutine reada(rekord,lancuch,wartosc,default)
320       implicit none
321       character*(*) rekord,lancuch
322       double precision wartosc,default
323       integer ilen,iread
324       external ilen
325       iread=index(rekord,lancuch)
326       if (iread.eq.0) then
327         wartosc=default 
328         return
329       endif   
330       iread=iread+ilen(lancuch)+1
331       read (rekord(iread:),*) wartosc
332       return
333       end
334 c----------------------------------------------------------------------------
335       subroutine multreada(rekord,lancuch,tablica,dim,default)
336       implicit none
337       integer dim,i
338       double precision tablica(dim),default
339       character*(*) rekord,lancuch
340       integer ilen,iread
341       external ilen
342       do i=1,dim
343         tablica(i)=default 
344       enddo
345       iread=index(rekord,lancuch)
346       if (iread.eq.0) return
347       iread=iread+ilen(lancuch)+1
348       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
349    10 return
350       end
351 c----------------------------------------------------------------------------
352       subroutine readi(rekord,lancuch,wartosc,default)
353       implicit none
354       character*(*) rekord,lancuch
355       integer wartosc,default
356       integer ilen,iread
357       external ilen
358       iread=index(rekord,lancuch)
359       if (iread.eq.0) then
360         wartosc=default 
361         return
362       endif   
363       iread=iread+ilen(lancuch)+1
364       read (rekord(iread:),*) wartosc
365       return
366       end
367 c----------------------------------------------------------------------------
368       subroutine split_string(rekord,tablica,dim,nsub)
369       implicit none
370       integer dim,nsub,i,ii,ll,kk
371       character*(*) tablica(dim)
372       character*(*) rekord
373       integer ilen
374       external ilen
375       do i=1,dim
376         tablica(i)=" "
377       enddo
378       ii=1
379       ll = ilen(rekord)
380       nsub=0
381       do i=1,dim
382 C Find the start of term name
383         kk = 0
384         do while (ii.le.ll .and. rekord(ii:ii).eq." ")
385           ii = ii+1
386         enddo
387 C Parse the name into TABLICA(i) until blank found
388         do while (ii.le.ll .and. rekord(ii:ii).ne." ")
389           kk = kk+1
390           tablica(i)(kk:kk)=rekord(ii:ii)
391           ii = ii+1
392         enddo
393         if (kk.gt.0) nsub=nsub+1
394         if (ii.gt.ll) return
395       enddo
396       return
397       end
398 c----------------------------------------------------------------------------
399       subroutine card_concat(card)
400       include 'DIMENSIONS'
401       include 'COMMON.IOUNITS'
402       character*(*) card
403       character*80 karta,ucase
404       external ilen
405       read (inp,'(a)') karta
406       karta=ucase(karta)
407       card=' '
408       do while (karta(80:80).eq.'&')
409         card=card(:ilen(card)+1)//karta(:79)
410         read (inp,'(a)') karta
411         karta=ucase(karta)
412       enddo
413       card=card(:ilen(card)+1)//karta
414       return
415       end
416 c----------------------------------------------------------------------------
417       subroutine openunits
418       include 'DIMENSIONS'    
419       include 'COMMON.IOUNITS'
420       include 'COMMON.CONTROL'
421       integer lenpre,lenpot,ilen
422       external ilen
423       character*16 cfrom_pdb,cprint
424       character*16 ucase
425       call getenv('INPUT',prefix)
426       call getenv('OUTPUT',prefout)
427       call getenv('INTIN',prefintin)
428       call getenv('PDB',cfrom_pdb)
429       call getenv('PRINTCOOR',cprint)
430       from_cart = index(ucase(cfrom_pdb),'CART').gt.0
431       from_cx = index(ucase(cfrom_pdb),'CX').gt.0
432       lprint_cart = index(ucase(cprint),'CART').gt.0
433       lprint_int = index(ucase(cprint),'INT').gt.0
434 c      if (.not.lprint_cart .and. .not.lprint_int) lprint_int=.true.
435       lenpre=ilen(prefix)
436       lenout=ilen(prefout)
437       lenint=ilen(prefintin)
438 C Get the names and open the input files
439       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
440       outname=prefout(:lenout)//'_clust.out'
441       if (from_cart) then
442         intinname=prefintin(:lenint)//'.x'
443       else if (from_cx) then
444         write (iout,*) "cx files: ",prefintin(:ilen(prefintin))
445         call split_string(prefintin,cxfiles(1),maxfiles,nfiles)
446         write (iout,*) "nfiles",nfiles
447         write (iout,*) "Split cxfiles"
448         do i=1,nfiles
449           cxfiles(i)=cxfiles(i)(:ilen(cxfiles(i)))//'.cx'
450           write (iout,*) cxfiles(i)(:ilen(cxfiles(i)))
451         enddo
452       else
453         intinname=prefintin(:lenint)//'.int'
454       endif
455       if (lprint_cart) then
456         if (.not. from_cx) then
457           intname=prefintin(:lenint)//'_clust'//'.x'
458         else
459           intname=cxfiles(1)(:ilen(cxfiles(1)))//'_clust'//'.x'
460         endif
461       else if (lprint_int) then
462         if (.not. from_cx) then
463           intname=prefintin(:lenint)//'_clust'//'.int'
464         else
465           intname=cxfiles(1)(:ilen(cxfiles(1)))//'_clust'//'.int'
466         endif
467       endif
468       rmsname=prefintin(:lenint)//'.rms'
469       statinname=prefintin(:lenint)//'.stat'
470       print *,statinname
471       statoutname=prefintin(:lenint)//'_clust.stat'
472       open (jplot,file=prefout(:ilen(prefout))//'.tex',
473      &       status='unknown')
474       print *,'unit',jplot,' opened'
475       if (.not. from_cx) open (intin,file=intinname,status='old')
476       print *,'unit',intin,' opened'
477       open (jrms,file=rmsname,status='unknown')
478       open (jstatin,file=statinname,status='unknown')
479       open (jstatout,file=statoutname,status='unknown')
480 #ifdef AIX
481       open(iout,file=outname,status='unknown')
482       print *,'unit',iout,' opened'
483       open(igeom,file=intname,status='unknown',position='append')
484       print *,'unit',igeom,' opened'
485 #elif defined(G77)
486       open(iout,file=outname,status='unknown')
487       open(igeom,file=intname,status='unknown',access='append')
488 #else
489       open(iout,file=outname,status='unknown')
490       open(igeom,file=intname,status='unknown',position='append')
491 #endif
492       return
493       end