xdrf added to cluster
[unres.git] / source / cluster / wham / src-M / readrtns.F
1       subroutine read_control
2 C
3 C Read molecular data
4 C
5       implicit none
6       include 'DIMENSIONS'
7       include 'sizesclu.dat'
8       include 'COMMON.IOUNITS'
9       include 'COMMON.TIME1'
10       include 'COMMON.SBRIDGE'
11       include 'COMMON.CONTROL'
12       include 'COMMON.CLUSTER'
13       include 'COMMON.CHAIN'
14       include 'COMMON.HEADER'
15       include 'COMMON.FFIELD'
16       include 'COMMON.FREE'
17       include 'COMMON.INTERACT'
18       character*320 controlcard,ucase
19 #ifdef MPL
20       include 'COMMON.INFO'
21 #endif
22       integer i
23
24       read (INP,'(a80)') titel
25       call card_concat(controlcard)
26
27       call readi(controlcard,'NRES',nres,0)
28       call readi(controlcard,'RESCALE',rescale_mode,2)
29       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
30       write (iout,*) "DISTCHAINMAX",distchainmax
31       call readi(controlcard,'PDBOUT',outpdb,0)
32       call readi(controlcard,'MOL2OUT',outmol2,0)
33       refstr=(index(controlcard,'REFSTR').gt.0)
34       write (iout,*) "REFSTR",refstr
35       pdbref=(index(controlcard,'PDBREF').gt.0)
36       iscode=index(controlcard,'ONE_LETTER')
37       tree=(index(controlcard,'MAKE_TREE').gt.0)
38       min_var=(index(controlcard,'MINVAR').gt.0)
39       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
40       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
41       call readi(controlcard,'NCUT',ncut,1)
42       call readi(controlcard,'SYM',symetr,1)
43       write (iout,*) 'sym', symetr
44       call readi(controlcard,'NSTART',nstart,0)
45       call readi(controlcard,'NEND',nend,0)
46       call reada(controlcard,'ECUT',ecut,10.0d0)
47       call reada(controlcard,'PROB',prob_limit,0.99d0)
48       write (iout,*) "Probability limit",prob_limit
49       lgrp=(index(controlcard,'LGRP').gt.0)
50       caonly=(index(controlcard,'CA_ONLY').gt.0)
51       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
52       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
53       call readi(controlcard,'IOPT',iopt,2) 
54       lside = index(controlcard,"SIDE").gt.0
55       efree = index(controlcard,"EFREE").gt.0
56       call readi(controlcard,'NTEMP',nT,1)
57       write (iout,*) "nT",nT
58       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
59       write (iout,*) "nT",nT
60       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
61       do i=1,nT
62         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
63       enddo
64       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
65       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
66       lprint_int=index(controlcard,"PRINT_INT") .gt.0
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       implicit none
76       include 'DIMENSIONS'
77       include 'COMMON.IOUNITS'
78       include 'COMMON.GEO'
79       include 'COMMON.VAR'
80       include 'COMMON.INTERACT'
81       include 'COMMON.LOCAL'
82       include 'COMMON.NAMES'
83       include 'COMMON.CHAIN'
84       include 'COMMON.FFIELD'
85       include 'COMMON.SBRIDGE'
86       include 'COMMON.HEADER'
87       include 'COMMON.CONTROL'
88       include 'COMMON.CONTACTS'
89       include 'COMMON.TIME1'
90 #ifdef MPL
91       include 'COMMON.INFO'
92 #endif
93       character*4 sequence(maxres)
94       character*800 weightcard
95       integer rescode
96       double precision x(maxvar)
97       integer itype_pdb(maxres)
98       logical seq_comp
99       integer i,j
100 C
101 C Body
102 C
103 C Read weights of the subsequent energy terms.
104       call card_concat(weightcard)
105       call reada(weightcard,'WSC',wsc,1.0d0)
106       call reada(weightcard,'WLONG',wsc,wsc)
107       call reada(weightcard,'WSCP',wscp,1.0d0)
108       call reada(weightcard,'WELEC',welec,1.0D0)
109       call reada(weightcard,'WVDWPP',wvdwpp,welec)
110       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
111       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
112       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
113       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
114       call reada(weightcard,'WTURN3',wturn3,1.0D0)
115       call reada(weightcard,'WTURN4',wturn4,1.0D0)
116       call reada(weightcard,'WTURN6',wturn6,1.0D0)
117       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
118       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
119       call reada(weightcard,'WBOND',wbond,1.0D0)
120       call reada(weightcard,'WTOR',wtor,1.0D0)
121       call reada(weightcard,'WTORD',wtor_d,1.0D0)
122       call reada(weightcard,'WANG',wang,1.0D0)
123       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
124       call reada(weightcard,'SCAL14',scal14,0.4D0)
125       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
126       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
127       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
128       if (index(weightcard,'SOFT').gt.0) ipot=6
129 C 12/1/95 Added weight for the multi-body term WCORR
130       call reada(weightcard,'WCORRH',wcorr,1.0D0)
131       if (wcorr4.gt.0.0d0) wcorr=wcorr4
132       weights(1)=wsc
133       weights(2)=wscp
134       weights(3)=welec
135       weights(4)=wcorr
136       weights(5)=wcorr5
137       weights(6)=wcorr6
138       weights(7)=wel_loc
139       weights(8)=wturn3
140       weights(9)=wturn4
141       weights(10)=wturn6
142       weights(11)=wang
143       weights(12)=wscloc
144       weights(13)=wtor
145       weights(14)=wtor_d
146       weights(15)=wstrain
147       weights(16)=wvdwpp
148       weights(17)=wbond
149       weights(18)=scal14
150       weights(19)=wsccor
151       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
152      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
153      &  wturn4,wturn6,wsccor
154    10 format (/'Energy-term weights (unscaled):'//
155      & 'WSCC=   ',f10.6,' (SC-SC)'/
156      & 'WSCP=   ',f10.6,' (SC-p)'/
157      & 'WELEC=  ',f10.6,' (p-p electr)'/
158      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
159      & 'WBOND=  ',f10.6,' (stretching)'/
160      & 'WANG=   ',f10.6,' (bending)'/
161      & 'WSCLOC= ',f10.6,' (SC local)'/
162      & 'WTOR=   ',f10.6,' (torsional)'/
163      & 'WTORD=  ',f10.6,' (double torsional)'/
164      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
165      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
166      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
167      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
168      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
169      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
170      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
171      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
172      & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
173
174       if (wcorr4.gt.0.0d0) then
175         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
176      &   'between contact pairs of peptide groups'
177         write (iout,'(2(a,f5.3/))')
178      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
179      &  'Range of quenching the correlation terms:',2*delt_corr
180       else if (wcorr.gt.0.0d0) then
181         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
182      &   'between contact pairs of peptide groups'
183       endif
184       write (iout,'(a,f8.3)')
185      &  'Scaling factor of 1,4 SC-p interactions:',scal14
186       write (iout,'(a,f8.3)')
187      &  'General scaling factor of SC-p interactions:',scalscp
188       r0_corr=cutoff_corr-delt_corr
189       do i=1,20
190         aad(i,1)=scalscp*aad(i,1)
191         aad(i,2)=scalscp*aad(i,2)
192         bad(i,1)=scalscp*bad(i,1)
193         bad(i,2)=scalscp*bad(i,2)
194       enddo
195
196       call flush(iout)
197       print *,'indpdb=',indpdb,' pdbref=',pdbref
198
199 C Read sequence if not taken from the pdb file.
200       if (iscode.gt.0) then
201         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
202       else
203         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
204       endif
205 C Convert sequence to numeric code
206       do i=1,nres
207         itype(i)=rescode(i,sequence(i),iscode)
208       enddo
209       print *,nres
210       print '(20i4)',(itype(i),i=1,nres)
211
212       do i=1,nres
213 #ifdef PROCOR
214         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
215 #else
216         if (itype(i).eq.21) then
217 #endif
218           itel(i)=0
219 #ifdef PROCOR
220         else if (itype(i+1).ne.20) then
221 #else
222         else if (itype(i).ne.20) then
223 #endif
224           itel(i)=1
225         else
226           itel(i)=2
227         endif
228       enddo
229       write (iout,*) "ITEL"
230       do i=1,nres-1
231         write (iout,*) i,itype(i),itel(i)
232       enddo
233
234       print *,'Call Read_Bridge.'
235       call read_bridge
236       nnt=1
237       nct=nres
238       print *,'NNT=',NNT,' NCT=',NCT
239       if (itype(1).eq.21) nnt=2
240       if (itype(nres).eq.21) nct=nct-1
241       if (nstart.lt.nnt) nstart=nnt
242       if (nend.gt.nct .or. nend.eq.0) nend=nct
243       write (iout,*) "nstart",nstart," nend",nend
244       nres0=nres
245 c      if (pdbref) then
246 c        read(inp,'(a)') pdbfile
247 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
248 c        open(ipdbin,file=pdbfile,status='old',err=33)
249 c        goto 34 
250 c  33    write (iout,'(a)') 'Error opening PDB file.'
251 c        stop
252 c  34    continue
253 c        print *,'Begin reading pdb data'
254 c        call readpdb
255 c        print *,'Finished reading pdb data'
256 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
257 c        do i=1,nres
258 c          itype_pdb(i)=itype(i)
259 c        enddo
260 c        close (ipdbin)
261 c        write (iout,'(a,i3)') 'nsup=',nsup
262 c        nstart_seq=nnt
263 c        if (nsup.le.(nct-nnt+1)) then
264 c          do i=0,nct-nnt+1-nsup
265 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
266 c              nstart_seq=nnt+i
267 c              goto 111
268 c            endif
269 c          enddo
270 c          write (iout,'(a)') 
271 c     &            'Error - sequences to be superposed do not match.'
272 c          stop
273 c        else
274 c          do i=0,nsup-(nct-nnt+1)
275 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) 
276 c     &      then
277 c              nstart_sup=nstart_sup+i
278 c              nsup=nct-nnt+1
279 c              goto 111
280 c            endif
281 c          enddo 
282 c          write (iout,'(a)') 
283 c     &            'Error - sequences to be superposed do not match.'
284 c        endif
285 c  111   continue
286 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
287 c     &                 ' nstart_seq=',nstart_seq
288 c      endif
289       call init_int_table
290       call setup_var
291       write (iout,*) "molread: REFSTR",refstr
292       if (refstr) then
293         if (.not.pdbref) then
294           call read_angles(inp,*38)
295           goto 39
296    38     write (iout,'(a)') 'Error reading reference structure.'
297 #ifdef MPL
298           call mp_stopall(Error_Msg)
299 #else
300           stop 'Error reading reference structure'
301 #endif
302    39     call chainbuild     
303           nstart_sup=nnt
304           nstart_seq=nnt
305           nsup=nct-nnt+1
306           do i=1,2*nres
307             do j=1,3
308               cref(j,i)=c(j,i)
309             enddo
310           enddo
311         endif
312         call contact(.true.,ncont_ref,icont_ref)
313       endif
314       return
315       end
316 c-----------------------------------------------------------------------------
317       logical function seq_comp(itypea,itypeb,length)
318       implicit none
319       integer length,itypea(length),itypeb(length)
320       integer i
321       do i=1,length
322         if (itypea(i).ne.itypeb(i)) then
323           seq_comp=.false.
324           return
325         endif
326       enddo
327       seq_comp=.true.
328       return
329       end
330 c-----------------------------------------------------------------------------
331       subroutine read_bridge
332 C Read information about disulfide bridges.
333       implicit none
334       include 'DIMENSIONS'
335       include 'COMMON.IOUNITS'
336       include 'COMMON.GEO'
337       include 'COMMON.VAR'
338       include 'COMMON.INTERACT'
339       include 'COMMON.LOCAL'
340       include 'COMMON.NAMES'
341       include 'COMMON.CHAIN'
342       include 'COMMON.FFIELD'
343       include 'COMMON.SBRIDGE'
344       include 'COMMON.HEADER'
345       include 'COMMON.CONTROL'
346       include 'COMMON.TIME1'
347 #ifdef MPL
348       include 'COMMON.INFO'
349 #endif
350       integer i,j
351 C Read bridging residues.
352       read (inp,*) ns,(iss(i),i=1,ns)
353       print *,'ns=',ns
354 C Check whether the specified bridging residues are cystines.
355       do i=1,ns
356         if (itype(iss(i)).ne.1) then
357           write (iout,'(2a,i3,a)') 
358      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
359      &   ' can form a disulfide bridge?!!!'
360           write (*,'(2a,i3,a)') 
361      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
362      &   ' can form a disulfide bridge?!!!'
363 #ifdef MPL
364          call mp_stopall(error_msg)
365 #else
366          stop
367 #endif
368         endif
369       enddo
370 C Read preformed bridges.
371       if (ns.gt.0) then
372       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
373       if (nss.gt.0) then
374         nhpb=nss
375 C Check if the residues involved in bridges are in the specified list of
376 C bridging residues.
377         do i=1,nss
378           do j=1,i-1
379             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
380      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
381               write (iout,'(a,i3,a)') 'Disulfide pair',i,
382      &      ' contains residues present in other pairs.'
383               write (*,'(a,i3,a)') 'Disulfide pair',i,
384      &      ' contains residues present in other pairs.'
385 #ifdef MPL
386               call mp_stopall(error_msg)
387 #else
388               stop 
389 #endif
390             endif
391           enddo
392           do j=1,ns
393             if (ihpb(i).eq.iss(j)) goto 10
394           enddo
395           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
396    10     continue
397           do j=1,ns
398             if (jhpb(i).eq.iss(j)) goto 20
399           enddo
400           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
401    20     continue
402           dhpb(i)=dbr
403           forcon(i)=fbr
404         enddo
405         do i=1,nss
406           ihpb(i)=ihpb(i)+nres
407           jhpb(i)=jhpb(i)+nres
408         enddo
409       endif
410       endif
411       return
412       end
413 c----------------------------------------------------------------------------
414       subroutine read_angles(kanal,*)
415       implicit none
416       include 'DIMENSIONS'
417       include 'COMMON.GEO'
418       include 'COMMON.VAR'
419       include 'COMMON.CHAIN'
420       include 'COMMON.IOUNITS'
421       integer i,kanal
422       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
423       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
424       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
425       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
426       do i=1,nres
427         theta(i)=deg2rad*theta(i)
428         phi(i)=deg2rad*phi(i)
429         alph(i)=deg2rad*alph(i)
430         omeg(i)=deg2rad*omeg(i)
431       enddo
432       return
433    10 return1
434       end
435 c----------------------------------------------------------------------------
436       subroutine reada(rekord,lancuch,wartosc,default)
437       implicit none
438       character*(*) rekord,lancuch
439       double precision wartosc,default
440       integer ilen,iread
441       external ilen
442       iread=index(rekord,lancuch)
443       if (iread.eq.0) then
444         wartosc=default 
445         return
446       endif   
447       iread=iread+ilen(lancuch)+1
448       read (rekord(iread:),*) wartosc
449       return
450       end
451 c----------------------------------------------------------------------------
452       subroutine multreada(rekord,lancuch,tablica,dim,default)
453       implicit none
454       integer dim,i
455       double precision tablica(dim),default
456       character*(*) rekord,lancuch
457       integer ilen,iread
458       external ilen
459       do i=1,dim
460         tablica(i)=default 
461       enddo
462       iread=index(rekord,lancuch)
463       if (iread.eq.0) return
464       iread=iread+ilen(lancuch)+1
465       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
466    10 return
467       end
468 c----------------------------------------------------------------------------
469       subroutine readi(rekord,lancuch,wartosc,default)
470       implicit none
471       character*(*) rekord,lancuch
472       integer wartosc,default
473       integer ilen,iread
474       external ilen
475       iread=index(rekord,lancuch)
476       if (iread.eq.0) then
477         wartosc=default 
478         return
479       endif   
480       iread=iread+ilen(lancuch)+1
481       read (rekord(iread:),*) wartosc
482       return
483       end
484 c----------------------------------------------------------------------------
485       subroutine card_concat(card)
486       include 'DIMENSIONS'
487       include 'COMMON.IOUNITS'
488       character*(*) card
489       character*80 karta,ucase
490       external ilen
491       read (inp,'(a)') karta
492       karta=ucase(karta)
493       card=' '
494       do while (karta(80:80).eq.'&')
495         card=card(:ilen(card)+1)//karta(:79)
496         read (inp,'(a)') karta
497         karta=ucase(karta)
498       enddo
499       card=card(:ilen(card)+1)//karta
500       return
501       end
502 c----------------------------------------------------------------------------
503       subroutine openunits
504       implicit none
505       include 'DIMENSIONS'    
506 #ifdef MPI
507       include "mpif.h"
508       character*3 liczba
509       include "COMMON.MPI"
510 #endif
511       include 'COMMON.IOUNITS'
512       include 'COMMON.CONTROL'
513       integer lenpre,lenpot,ilen
514       external ilen
515       character*16 cformat,cprint
516       character*16 ucase
517       integer lenint,lenout
518       call getenv('INPUT',prefix)
519       call getenv('OUTPUT',prefout)
520       call getenv('INTIN',prefintin)
521       call getenv('COORD',cformat)
522       call getenv('PRINTCOOR',cprint)
523       call getenv('SCRATCHDIR',scratchdir)
524       from_bx=.true.
525       from_cx=.false.
526       if (index(ucase(cformat),'CX').gt.0) then
527         from_cx=.true.
528         from_bx=.false.
529       endif
530       from_cart=.true.
531       lenpre=ilen(prefix)
532       lenout=ilen(prefout)
533       lenint=ilen(prefintin)
534 C Get the names and open the input files
535       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
536 #ifdef MPI
537       write (liczba,'(bz,i3.3)') me
538       outname=prefout(:lenout)//'_clust.out_'//liczba
539 #else
540       outname=prefout(:lenout)//'_clust.out'
541 #endif
542       if (from_bx) then
543         intinname=prefintin(:lenint)//'.bx'
544       else if (from_cx) then
545         intinname=prefintin(:lenint)//'.cx'
546       else
547         intinname=prefintin(:lenint)//'.int'
548       endif
549       rmsname=prefintin(:lenint)//'.rms'
550       open (jplot,file=prefout(:ilen(prefout))//'.tex',
551      &       status='unknown')
552       open (jrms,file=rmsname,status='unknown')
553       open(iout,file=outname,status='unknown')
554 C Get parameter filenames and open the parameter files.
555       call getenv('BONDPAR',bondname)
556       open (ibond,file=bondname,status='old')
557       call getenv('THETPAR',thetname)
558       open (ithep,file=thetname,status='old')
559       call getenv('ROTPAR',rotname)
560       open (irotam,file=rotname,status='old')
561       call getenv('TORPAR',torname)
562       open (itorp,file=torname,status='old')
563       call getenv('TORDPAR',tordname)
564       open (itordp,file=tordname,status='old')
565       call getenv('FOURIER',fouriername)
566       open (ifourier,file=fouriername,status='old')
567       call getenv('ELEPAR',elename)
568       open (ielep,file=elename,status='old')
569       call getenv('SIDEPAR',sidename)
570       open (isidep,file=sidename,status='old')
571       call getenv('SIDEP',sidepname)
572       open (isidep1,file=sidepname,status="old")
573       call getenv('SCCORPAR',sccorname)
574       open (isccor,file=sccorname,status="old")
575 #ifndef OLDSCP
576 C
577 C 8/9/01 In the newest version SCp interaction constants are read from a file
578 C Use -DOLDSCP to use hard-coded constants instead.
579 C
580       call getenv('SCPPAR',scpname)
581       open (iscpp,file=scpname,status='old')
582 #endif
583       return
584       end