Rozgrzebany DIL
[unres.git] / source / unres / src_MD / parmread.F
1       subroutine parmread
2 C
3 C Read the parameters of the probability distributions of the virtual-bond
4 C valence angles and the side chains and energy parameters.
5 C
6 C Important! Energy-term weights ARE NOT read here; they are read from the
7 C main input file instead, because NO defaults have yet been set for these
8 C parameters.
9 C
10       implicit real*8 (a-h,o-z)
11       include 'DIMENSIONS'
12 #ifdef MPI
13       include "mpif.h"
14       integer IERROR
15 #endif
16       include 'COMMON.IOUNITS'
17       include 'COMMON.CHAIN'
18       include 'COMMON.INTERACT'
19       include 'COMMON.GEO'
20       include 'COMMON.LOCAL'
21       include 'COMMON.TORSION'
22       include 'COMMON.SCCOR'
23       include 'COMMON.SCROT'
24       include 'COMMON.FFIELD'
25       include 'COMMON.NAMES'
26       include 'COMMON.SBRIDGE'
27       include 'COMMON.MD'
28       include 'COMMON.SETUP'
29       character*1 t1,t2,t3
30       character*1 onelett(4) /"G","A","P","D"/
31       logical lprint,LaTeX
32       dimension blower(3,3,maxlob)
33       dimension b(13)
34       character*3 lancuch,ucase
35 C
36 C For printing parameters after they are read set the following in the UNRES
37 C C-shell script:
38 C
39 C setenv PRINT_PARM YES
40 C
41 C To print parameters in LaTeX format rather than as ASCII tables:
42 C
43 C setenv LATEX YES
44 C
45       call getenv_loc("PRINT_PARM",lancuch)
46       lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
47       call getenv_loc("LATEX",lancuch)
48       LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
49 C
50       dwa16=2.0d0**(1.0d0/6.0d0)
51       itypro=20
52 C Assign virtual-bond length
53       vbl=3.8D0
54       vblinv=1.0D0/vbl
55       vblinv2=vblinv*vblinv
56 c
57 c Read the virtual-bond parameters, masses, and moments of inertia
58 c and Stokes' radii of the peptide group and side chains
59 c
60 #ifdef CRYST_BOND
61       read (ibond,*) vbldp0,akp,mp,ip,pstok
62       do i=1,ntyp
63         nbondterm(i)=1
64         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
65         dsc(i) = vbldsc0(1,i)
66         if (i.eq.10) then
67           dsc_inv(i)=0.0D0
68         else
69           dsc_inv(i)=1.0D0/dsc(i)
70         endif
71       enddo
72 #else
73       read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
74       do i=1,ntyp
75         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
76      &   j=1,nbondterm(i)),msc(i),isc(i),restok(i)
77         dsc(i) = vbldsc0(1,i)
78         if (i.eq.10) then
79           dsc_inv(i)=0.0D0
80         else
81           dsc_inv(i)=1.0D0/dsc(i)
82         endif
83       enddo
84 #endif
85       if (lprint) then
86         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
87         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
88      &   'inertia','Pstok'
89         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
90         do i=1,ntyp
91           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
92      &      vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
93           do j=2,nbondterm(i)
94             write (iout,'(13x,3f10.5)')
95      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
96           enddo
97         enddo
98       endif
99 #ifdef CRYST_THETA
100 C
101 C Read the parameters of the probability distribution/energy expression 
102 C of the virtual-bond valence angles theta
103 C
104       do i=1,ntyp
105         read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
106      &    (bthet(j,i),j=1,2)
107         read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
108         read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
109         read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
110         sigc0(i)=sigc0(i)**2
111       enddo
112       close (ithep)
113       if (lprint) then
114       if (.not.LaTeX) then
115         write (iout,'(a)') 
116      &    'Parameters of the virtual-bond valence angles:'
117         write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
118      & '    ATHETA0   ','         A1   ','        A2    ',
119      & '        B1    ','         B2   '        
120         do i=1,ntyp
121           write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
122      &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
123         enddo
124         write (iout,'(/a/9x,5a/79(1h-))') 
125      & 'Parameters of the expression for sigma(theta_c):',
126      & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
127      & '     ALPH3    ','    SIGMA0C   '        
128         do i=1,ntyp
129           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
130      &      (polthet(j,i),j=0,3),sigc0(i) 
131         enddo
132         write (iout,'(/a/9x,5a/79(1h-))') 
133      & 'Parameters of the second gaussian:',
134      & '    THETA0    ','     SIGMA0   ','        G1    ',
135      & '        G2    ','         G3   '        
136         do i=1,ntyp
137           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
138      &       sig0(i),(gthet(j,i),j=1,3)
139         enddo
140        else
141         write (iout,'(a)') 
142      &    'Parameters of the virtual-bond valence angles:'
143         write (iout,'(/a/9x,5a/79(1h-))') 
144      & 'Coefficients of expansion',
145      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
146      & '   b1*10^1    ','    b2*10^1   '        
147         do i=1,ntyp
148           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
149      &        a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
150         enddo
151         write (iout,'(/a/9x,5a/79(1h-))') 
152      & 'Parameters of the expression for sigma(theta_c):',
153      & ' alpha0       ','  alph1       ',' alph2        ',
154      & ' alhp3        ','   sigma0c    '        
155         do i=1,ntyp
156           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
157      &      (polthet(j,i),j=0,3),sigc0(i) 
158         enddo
159         write (iout,'(/a/9x,5a/79(1h-))') 
160      & 'Parameters of the second gaussian:',
161      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
162      & '        G2    ','   G3*10^1    '        
163         do i=1,ntyp
164           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
165      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
166         enddo
167       endif
168       endif
169 #else 
170
171 C Read the parameters of Utheta determined from ab initio surfaces
172 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
173 C
174       read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
175      &  ntheterm3,nsingle,ndouble
176       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
177       read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
178       do i=1,maxthetyp
179         do j=1,maxthetyp
180           do k=1,maxthetyp
181             aa0thet(i,j,k)=0.0d0
182             do l=1,ntheterm
183               aathet(l,i,j,k)=0.0d0
184             enddo
185             do l=1,ntheterm2
186               do m=1,nsingle
187                 bbthet(m,l,i,j,k)=0.0d0
188                 ccthet(m,l,i,j,k)=0.0d0
189                 ddthet(m,l,i,j,k)=0.0d0
190                 eethet(m,l,i,j,k)=0.0d0
191               enddo
192             enddo
193             do l=1,ntheterm3
194               do m=1,ndouble
195                 do mm=1,ndouble
196                  ffthet(mm,m,l,i,j,k)=0.0d0
197                  ggthet(mm,m,l,i,j,k)=0.0d0
198                 enddo
199               enddo
200             enddo
201           enddo
202         enddo
203       enddo 
204       do i=1,nthetyp
205         do j=1,nthetyp
206           do k=1,nthetyp
207             read (ithep,'(3a)',end=111,err=111) res1,res2,res3
208             read (ithep,*,end=111,err=111) aa0thet(i,j,k)
209             read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
210             read (ithep,*,end=111,err=111)
211      &       ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
212      &        (ccthet(lll,ll,i,j,k),lll=1,nsingle),
213      &        (ddthet(lll,ll,i,j,k),lll=1,nsingle),
214      &        (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
215             read (ithep,*,end=111,err=111)
216      &      (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
217      &         ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
218      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
219           enddo
220         enddo
221       enddo
222 C
223 C For dummy ends assign glycine-type coefficients of theta-only terms; the
224 C coefficients of theta-and-gamma-dependent terms are zero.
225 C
226       do i=1,nthetyp
227         do j=1,nthetyp
228           do l=1,ntheterm
229             aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
230             aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
231           enddo
232           aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
233           aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
234         enddo
235         do l=1,ntheterm
236           aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
237         enddo
238         aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
239       enddo
240 C
241 C Control printout of the coefficients of virtual-bond-angle potentials
242 C
243       if (lprint) then
244         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
245         do i=1,nthetyp+1
246           do j=1,nthetyp+1
247             do k=1,nthetyp+1
248               write (iout,'(//4a)') 
249      &         'Type ',onelett(i),onelett(j),onelett(k) 
250               write (iout,'(//a,10x,a)') " l","a[l]"
251               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
252               write (iout,'(i2,1pe15.5)')
253      &           (l,aathet(l,i,j,k),l=1,ntheterm)
254             do l=1,ntheterm2
255               write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') 
256      &          "b",l,"c",l,"d",l,"e",l
257               do m=1,nsingle
258                 write (iout,'(i2,4(1pe15.5))') m,
259      &          bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
260      &          ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
261               enddo
262             enddo
263             do l=1,ntheterm3
264               write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
265      &          "f+",l,"f-",l,"g+",l,"g-",l
266               do m=2,ndouble
267                 do n=1,m-1
268                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
269      &              ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
270      &              ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
271                 enddo
272               enddo
273             enddo
274           enddo
275         enddo
276       enddo
277       call flush(iout)
278       endif
279       write (2,*) "Start reading THETA_PDB"
280       do i=1,ntyp
281         read (ithep_pdb,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
282      &    (bthet(j,i),j=1,2)
283         read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
284         read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
285         read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
286         sigc0(i)=sigc0(i)**2
287       enddo
288       write (2,*) "End reading THETA_PDB"
289       close (ithep_pdb)
290 #endif
291       close(ithep)
292 #ifdef CRYST_SC
293 C
294 C Read the parameters of the probability distribution/energy expression
295 C of the side chains.
296 C
297       do i=1,ntyp
298         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
299         if (i.eq.10) then
300           dsc_inv(i)=0.0D0
301         else
302           dsc_inv(i)=1.0D0/dsc(i)
303         endif
304         if (i.ne.10) then
305         do j=1,nlob(i)
306           do k=1,3
307             do l=1,3
308               blower(l,k,j)=0.0D0
309             enddo
310           enddo
311         enddo  
312         bsc(1,i)=0.0D0
313         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
314      &    ((blower(k,l,1),l=1,k),k=1,3)
315         censc(1,1,-i)=censc(1,1,i)
316         censc(2,1,-i)=censc(2,1,i)
317         censc(3,1,-i)=-censc(3,1,i)
318         
319         do j=2,nlob(i)
320           read (irotam,*,end=112,err=112) bsc(j,i)
321           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
322      &                                 ((blower(k,l,j),l=1,k),k=1,3)
323         censc(1,j,-i)=censc(1,j,i)
324         censc(2,j,-i)=censc(2,j,i)
325         censc(3,j,-i)=-censc(3,j,i)
326 C BSC is amplitude of Gaussian
327         enddo
328         do j=1,nlob(i)
329           do k=1,3
330             do l=1,k
331               akl=0.0D0
332               do m=1,3
333                 akl=akl+blower(k,m,j)*blower(l,m,j)
334               enddo
335               gaussc(k,l,j,i)=akl
336               gaussc(l,k,j,i)=akl
337               if (((k.eq.3).and.(l.ne.3))
338      &        .or.((l.eq.3).and.(k.ne.3))) then
339                 gaussc(k,l,j,-i)=-akl
340                 gaussc(l,k,j,-i)=-akl
341               else
342                 gaussc(k,l,j,-i)=akl
343                 gaussc(l,k,j,-i)=akl
344               endif
345             enddo
346           enddo 
347         enddo
348         endif
349       enddo
350       close (irotam)
351       if (lprint) then
352         write (iout,'(/a)') 'Parameters of side-chain local geometry'
353         do i=1,ntyp
354           nlobi=nlob(i)
355           if (nlobi.gt.0) then
356             if (LaTeX) then
357               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
358      &         ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
359                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
360      &                             'log h',(bsc(j,i),j=1,nlobi)
361                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
362      &        'x',((censc(k,j,i),k=1,3),j=1,nlobi)
363               do k=1,3
364                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
365      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
366               enddo
367             else
368               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
369               write (iout,'(a,f10.4,4(16x,f10.4))')
370      &                             'Center  ',(bsc(j,i),j=1,nlobi)
371               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
372      &           j=1,nlobi)
373               write (iout,'(a)')
374             endif
375           endif
376         enddo
377       endif
378 #else
379
380 C Read scrot parameters for potentials determined from all-atom AM1 calculations
381 C added by Urszula Kozlowska 07/11/2007
382 C
383       do i=1,ntyp
384         read (irotam,*,end=112,err=112) 
385        if (i.eq.10) then 
386          read (irotam,*,end=112,err=112) 
387        else
388          do j=1,65
389            read(irotam,*,end=112,err=112) sc_parmin(j,i)
390          enddo  
391        endif
392       enddo
393 C
394 C Read the parameters of the probability distribution/energy expression
395 C of the side chains.
396 C
397       do i=1,ntyp
398         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
399         if (i.eq.10) then
400           dsc_inv(i)=0.0D0
401         else
402           dsc_inv(i)=1.0D0/dsc(i)
403         endif
404         if (i.ne.10) then
405         do j=1,nlob(i)
406           do k=1,3
407             do l=1,3
408               blower(l,k,j)=0.0D0
409             enddo
410           enddo
411         enddo  
412         bsc(1,i)=0.0D0
413         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
414      &    ((blower(k,l,1),l=1,k),k=1,3)
415         do j=2,nlob(i)
416           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
417           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
418      &                                 ((blower(k,l,j),l=1,k),k=1,3)
419         enddo
420         do j=1,nlob(i)
421           do k=1,3
422             do l=1,k
423               akl=0.0D0
424               do m=1,3
425                 akl=akl+blower(k,m,j)*blower(l,m,j)
426               enddo
427               gaussc(k,l,j,i)=akl
428               gaussc(l,k,j,i)=akl
429             enddo
430           enddo 
431         enddo
432         endif
433       enddo
434       close (irotam_pdb)
435 #endif
436       close(irotam)
437
438 #ifdef CRYST_TOR
439 C
440 C Read torsional parameters in old format
441 C
442       read (itorp,*,end=113,err=113) ntortyp,nterm_old
443       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
444       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
445       do i=1,ntortyp
446         do j=1,ntortyp
447           read (itorp,'(a)')
448           do k=1,nterm_old
449             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i) 
450           enddo
451         enddo
452       enddo
453       close (itorp)
454       if (lprint) then
455         write (iout,'(/a/)') 'Torsional constants:'
456         do i=1,ntortyp
457           do j=1,ntortyp
458             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
459             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
460           enddo
461         enddo
462       endif
463 #else
464 C
465 C Read torsional parameters
466 C
467       read (itorp,*,end=113,err=113) ntortyp
468       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
469       do iblock=1,2
470       do i=-ntyp,-1
471        itortyp(i)=-itortyp(-i)
472       enddo
473 c      write (iout,*) 'ntortyp',ntortyp
474       do i=0,ntortyp-1
475         do j=-ntortyp,ntortyp
476           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
477      &          nlor(i,j,iblock)
478           nterm(-i,-j,iblock)=nterm(i,j,iblock)
479           nlor(-i,-j,iblock)=nlor(i,j,iblock)
480           v0ij=0.0d0
481           si=-1.0d0
482           do k=1,nterm(i,j,iblock)
483             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
484      &      v2(k,i,j,iblock) 
485             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
486             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
487             v0ij=v0ij+si*v1(k,i,j)
488             si=-si
489           enddo
490           do k=1,nlor(i,j,iblock)
491             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
492      &        vlor2(k,i,j),vlor3(k,i,j) 
493             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
494           enddo
495           v0(i,j,iblock)=v0ij
496           v0(-i,-j,iblock)=v0ij
497         enddo
498       enddo
499       enddo
500       close (itorp)
501       if (lprint) then
502         write (iout,'(/a/)') 'Torsional constants:'
503         do i=1,ntortyp
504           do j=1,ntortyp
505             write (iout,*) 'ityp',i,' jtyp',j
506             write (iout,*) 'Fourier constants'
507             do k=1,nterm(i,j)
508               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
509             enddo
510             write (iout,*) 'Lorenz constants'
511             do k=1,nlor(i,j)
512               write (iout,'(3(1pe15.5))') 
513      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
514             enddo
515           enddo
516         enddo
517       endif
518 C
519 C 6/23/01 Read parameters for double torsionals
520 C
521       do i=1,ntortyp
522         do j=1,ntortyp
523           do k=1,ntortyp
524             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
525             if (t1.ne.onelett(i) .or. t2.ne.onelett(j) 
526      &        .or. t3.ne.onelett(k)) then
527               write (iout,*) "Error in double torsional parameter file",
528      &         i,j,k,t1,t2,t3
529 #ifdef MPI
530               call MPI_Finalize(Ierror)
531 #endif
532                stop "Error in double torsional parameter file"
533             endif
534             read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
535      &         ntermd_2(i,j,k)
536             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
537      &         ntermd_1(i,j,k))
538             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
539      &         ntermd_1(i,j,k))
540             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
541      &         ntermd_1(i,j,k))
542             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
543      &         ntermd_1(i,j,k))
544             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
545      &         v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
546      &         m=1,l-1),l=1,ntermd_2(i,j,k))
547           enddo
548         enddo
549       enddo
550       if (lprint) then
551       write (iout,*) 
552       write (iout,*) 'Constants for double torsionals'
553       do i=1,ntortyp
554         do j=1,ntortyp 
555           do k=1,ntortyp
556             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
557      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
558             write (iout,*)
559             write (iout,*) 'Single angles:'
560             do l=1,ntermd_1(i,j,k)
561               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
562      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
563      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
564             enddo
565             write (iout,*)
566             write (iout,*) 'Pairs of angles:'
567             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
568             do l=1,ntermd_2(i,j,k)
569               write (iout,'(i5,20f10.5)') 
570      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
571             enddo
572             write (iout,*)
573             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
574             do l=1,ntermd_2(i,j,k)
575               write (iout,'(i5,20f10.5)') 
576      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
577             enddo
578             write (iout,*)
579           enddo
580         enddo
581       enddo
582       endif
583 #endif
584 C Read of Side-chain backbone correlation parameters
585 C Modified 11 May 2012 by Adasko
586 CCC
587 C
588       read (isccor,*,end=113,err=113) nsccortyp
589 #ifdef SCCORPDB
590       read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
591 #else
592       read (isccor,*,end=113,err=113) (isccortyp(i),i=-ntyp,ntyp)
593 #endif
594 c      write (iout,*) 'ntortyp',ntortyp
595       maxinter=3
596 cc maxinter is maximum interaction sites
597       do l=1,maxinter    
598       do i=1,nsccortyp
599         do j=1,nsccortyp
600           read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j)
601           v0ijsccor=0.0d0
602           si=-1.0d0
603   
604           do k=1,nterm_sccor(i,j)
605             read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
606      &    ,v2sccor(k,l,i,j) 
607             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
608             si=-si
609           enddo
610           do k=1,nlor_sccor(i,j)
611             read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
612      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j) 
613             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
614      &(1+vlor3sccor(k,i,j)**2)
615           enddo
616           v0sccor(i,j)=v0ijsccor
617         enddo
618       enddo
619       enddo
620       close (isccor)
621       
622       if (lprint) then
623         write (iout,'(/a/)') 'Torsional constants:'
624         do i=1,nsccortyp
625           do j=1,nsccortyp
626             write (iout,*) 'ityp',i,' jtyp',j
627             write (iout,*) 'Fourier constants'
628             do k=1,nterm_sccor(i,j)
629       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
630             enddo
631             write (iout,*) 'Lorenz constants'
632             do k=1,nlor_sccor(i,j)
633               write (iout,'(3(1pe15.5))') 
634      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
635             enddo
636           enddo
637         enddo
638       endif
639 C
640 C
641 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
642 C         interaction energy of the Gly, Ala, and Pro prototypes.
643 C
644       if (lprint) then
645         write (iout,*)
646         write (iout,*) "Coefficients of the cumulants"
647       endif
648       read (ifourier,*) nloctyp
649       do i=1,nloctyp
650         read (ifourier,*,end=115,err=115)
651         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
652         if (lprint) then
653         write (iout,*) 'Type',i
654         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
655         endif
656         B1(1,i)  = b(3)
657         B1(2,i)  = b(5)
658         B1(1,-i) = b(3)
659         B1(2,-i) = -b(5)
660 c        b1(1,i)=0.0d0
661 c        b1(2,i)=0.0d0
662         B1tilde(1,i) = b(3)
663         B1tilde(2,i) =-b(5)
664         B1tilde(1,-i) =-b(3)
665         B1tilde(2,-i) =b(5) 
666 c        b1tilde(1,i)=0.0d0
667 c        b1tilde(2,i)=0.0d0
668         B2(1,i)  = b(2)
669         B2(2,i)  = b(4)
670         B2(1,-i)  =b(2)
671         B2(2,-i)  =-b(4)
672
673 c        b2(1,i)=0.0d0
674 c        b2(2,i)=0.0d0
675         CC(1,1,i)= b(7)
676         CC(2,2,i)=-b(7)
677         CC(2,1,i)= b(9)
678         CC(1,2,i)= b(9)
679         CC(1,1,-i)= b(7)
680         CC(2,2,-i)=-b(7)
681         CC(2,1,-i)=-b(9)
682         CC(1,2,-i)=-b(9)
683 c        CC(1,1,i)=0.0d0
684 c        CC(2,2,i)=0.0d0
685 c        CC(2,1,i)=0.0d0
686 c        CC(1,2,i)=0.0d0
687         Ctilde(1,1,i)=b(7)
688         Ctilde(1,2,i)=b(9)
689         Ctilde(2,1,i)=-b(9)
690         Ctilde(2,2,i)=b(7)
691         Ctilde(1,1,-i)=b(7)
692         Ctilde(1,2,-i)=-b(9)
693         Ctilde(2,1,-i)=b(9)
694         Ctilde(2,2,-i)=b(7)
695
696 c        Ctilde(1,1,i)=0.0d0
697 c        Ctilde(1,2,i)=0.0d0
698 c        Ctilde(2,1,i)=0.0d0
699 c        Ctilde(2,2,i)=0.0d0
700         DD(1,1,i)= b(6)
701         DD(2,2,i)=-b(6)
702         DD(2,1,i)= b(8)
703         DD(1,2,i)= b(8)
704         DD(1,1,-i)= b(6)
705         DD(2,2,-i)=-b(6)
706         DD(2,1,-i)=-b(8)
707         DD(1,2,-i)=-b(8)
708 c        DD(1,1,i)=0.0d0
709 c        DD(2,2,i)=0.0d0
710 c        DD(2,1,i)=0.0d0
711 c        DD(1,2,i)=0.0d0
712         Dtilde(1,1,i)=b(6)
713         Dtilde(1,2,i)=b(8)
714         Dtilde(2,1,i)=-b(8)
715         Dtilde(2,2,i)=b(6)
716         Dtilde(1,1,-i)=b(6)
717         Dtilde(1,2,-i)=-b(8)
718         Dtilde(2,1,-i)=b(8)
719         Dtilde(2,2,-i)=b(6)
720
721 c        Dtilde(1,1,i)=0.0d0
722 c        Dtilde(1,2,i)=0.0d0
723 c        Dtilde(2,1,i)=0.0d0
724 c        Dtilde(2,2,i)=0.0d0
725         EE(1,1,i)= b(10)+b(11)
726         EE(2,2,i)=-b(10)+b(11)
727         EE(2,1,i)= b(12)-b(13)
728         EE(1,2,i)= b(12)+b(13)
729         EE(1,1,-i)= b(10)+b(11)
730         EE(2,2,-i)=-b(10)+b(11)
731         EE(2,1,-i)=-b(12)+b(13)
732         EE(1,2,-i)=-b(12)-b(13)
733
734 c        ee(1,1,i)=1.0d0
735 c        ee(2,2,i)=1.0d0
736 c        ee(2,1,i)=0.0d0
737 c        ee(1,2,i)=0.0d0
738 c        ee(2,1,i)=ee(1,2,i)
739       enddo
740       if (lprint) then
741       do i=1,nloctyp
742         write (iout,*) 'Type',i
743         write (iout,*) 'B1'
744         write(iout,*) B1(1,i),B1(2,i)
745         write (iout,*) 'B2'
746         write(iout,*) B2(1,i),B2(2,i)
747         write (iout,*) 'CC'
748         do j=1,2
749           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
750         enddo
751         write(iout,*) 'DD'
752         do j=1,2
753           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
754         enddo
755         write(iout,*) 'EE'
756         do j=1,2
757           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
758         enddo
759       enddo
760       endif
761
762 C Read electrostatic-interaction parameters
763 C
764       if (lprint) then
765         write (iout,*)
766         write (iout,'(/a)') 'Electrostatic interaction constants:'
767         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
768      &            'IT','JT','APP','BPP','AEL6','AEL3'
769       endif
770       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
771       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
772       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
773       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
774       close (ielep)
775       do i=1,2
776         do j=1,2
777         rri=rpp(i,j)**6
778         app (i,j)=epp(i,j)*rri*rri 
779         bpp (i,j)=-2.0D0*epp(i,j)*rri
780         ael6(i,j)=elpp6(i,j)*4.2D0**6
781         ael3(i,j)=elpp3(i,j)*4.2D0**3
782         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
783      &                    ael6(i,j),ael3(i,j)
784         enddo
785       enddo
786 C
787 C Read side-chain interaction parameters.
788 C
789       read (isidep,*,end=117,err=117) ipot,expon
790       if (ipot.lt.1 .or. ipot.gt.5) then
791         write (iout,'(2a)') 'Error while reading SC interaction',
792      &               'potential file - unknown potential type.'
793 #ifdef MPI
794         call MPI_Finalize(Ierror)
795 #endif
796         stop
797       endif
798       expon2=expon/2
799       if(me.eq.king)
800      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
801      & ', exponents are ',expon,2*expon 
802       goto (10,20,30,30,40) ipot
803 C----------------------- LJ potential ---------------------------------
804    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
805      &   (sigma0(i),i=1,ntyp)
806       if (lprint) then
807         write (iout,'(/a/)') 'Parameters of the LJ potential:'
808         write (iout,'(a/)') 'The epsilon array:'
809         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
810         write (iout,'(/a)') 'One-body parameters:'
811         write (iout,'(a,4x,a)') 'residue','sigma'
812         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
813       endif
814       goto 50
815 C----------------------- LJK potential --------------------------------
816    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
817      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
818       if (lprint) then
819         write (iout,'(/a/)') 'Parameters of the LJK potential:'
820         write (iout,'(a/)') 'The epsilon array:'
821         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
822         write (iout,'(/a)') 'One-body parameters:'
823         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
824         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
825      &        i=1,ntyp)
826       endif
827       goto 50
828 C---------------------- GB or BP potential -----------------------------
829    30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
830      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
831      &  (alp(i),i=1,ntyp)
832 C For the GB potential convert sigma'**2 into chi'
833       if (ipot.eq.4) then
834         do i=1,ntyp
835           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
836         enddo
837       endif
838       if (lprint) then
839         write (iout,'(/a/)') 'Parameters of the BP potential:'
840         write (iout,'(a/)') 'The epsilon array:'
841         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
842         write (iout,'(/a)') 'One-body parameters:'
843         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
844      &       '    chip  ','    alph  '
845         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
846      &                     chip(i),alp(i),i=1,ntyp)
847       endif
848       goto 50
849 C--------------------- GBV potential -----------------------------------
850    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
851      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
852      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
853       if (lprint) then
854         write (iout,'(/a/)') 'Parameters of the GBV potential:'
855         write (iout,'(a/)') 'The epsilon array:'
856         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
857         write (iout,'(/a)') 'One-body parameters:'
858         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
859      &      's||/s_|_^2','    chip  ','    alph  '
860         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
861      &           sigii(i),chip(i),alp(i),i=1,ntyp)
862       endif
863    50 continue
864       close (isidep)
865 C-----------------------------------------------------------------------
866 C Calculate the "working" parameters of SC interactions.
867       do i=2,ntyp
868         do j=1,i-1
869           eps(i,j)=eps(j,i)
870         enddo
871       enddo
872       do i=1,ntyp
873         do j=i,ntyp
874           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
875           sigma(j,i)=sigma(i,j)
876           rs0(i,j)=dwa16*sigma(i,j)
877           rs0(j,i)=rs0(i,j)
878         enddo
879       enddo
880       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
881      & 'Working parameters of the SC interactions:',
882      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
883      & '  chi1   ','   chi2   ' 
884       do i=1,ntyp
885         do j=i,ntyp
886           epsij=eps(i,j)
887           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
888             rrij=sigma(i,j)
889           else
890             rrij=rr0(i)+rr0(j)
891           endif
892           r0(i,j)=rrij
893           r0(j,i)=rrij
894           rrij=rrij**expon
895           epsij=eps(i,j)
896           sigeps=dsign(1.0D0,epsij)
897           epsij=dabs(epsij)
898           aa(i,j)=epsij*rrij*rrij
899           bb(i,j)=-sigeps*epsij*rrij
900           aa(j,i)=aa(i,j)
901           bb(j,i)=bb(i,j)
902           if (ipot.gt.2) then
903             sigt1sq=sigma0(i)**2
904             sigt2sq=sigma0(j)**2
905             sigii1=sigii(i)
906             sigii2=sigii(j)
907             ratsig1=sigt2sq/sigt1sq
908             ratsig2=1.0D0/ratsig1
909             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
910             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
911             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
912           else
913             rsum_max=sigma(i,j)
914           endif
915 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
916             sigmaii(i,j)=rsum_max
917             sigmaii(j,i)=rsum_max 
918 c         else
919 c           sigmaii(i,j)=r0(i,j)
920 c           sigmaii(j,i)=r0(i,j)
921 c         endif
922 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
923           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
924             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
925             augm(i,j)=epsij*r_augm**(2*expon)
926 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
927             augm(j,i)=augm(i,j)
928           else
929             augm(i,j)=0.0D0
930             augm(j,i)=0.0D0
931           endif
932           if (lprint) then
933             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
934      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
935      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
936           endif
937         enddo
938       enddo
939 #ifdef OLDSCP
940 C
941 C Define the SC-p interaction constants (hard-coded; old style)
942 C
943       do i=1,20
944 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
945 C helix formation)
946 c       aad(i,1)=0.3D0*4.0D0**12
947 C Following line for constants currently implemented
948 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
949         aad(i,1)=1.5D0*4.0D0**12
950 c       aad(i,1)=0.17D0*5.6D0**12
951         aad(i,2)=aad(i,1)
952 C "Soft" SC-p repulsion
953         bad(i,1)=0.0D0
954 C Following line for constants currently implemented
955 c       aad(i,1)=0.3D0*4.0D0**6
956 C "Hard" SC-p repulsion
957         bad(i,1)=3.0D0*4.0D0**6
958 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
959         bad(i,2)=bad(i,1)
960 c       aad(i,1)=0.0D0
961 c       aad(i,2)=0.0D0
962 c       bad(i,1)=1228.8D0
963 c       bad(i,2)=1228.8D0
964       enddo
965 #else
966 C
967 C 8/9/01 Read the SC-p interaction constants from file
968 C
969       do i=1,ntyp
970         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
971       enddo
972       do i=1,ntyp
973         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
974         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
975         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
976         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
977       enddo
978
979       if (lprint) then
980         write (iout,*) "Parameters of SC-p interactions:"
981         do i=1,20
982           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
983      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
984         enddo
985       endif
986 #endif
987 C
988 C Define the constants of the disulfide bridge
989 C
990       ebr=-5.50D0
991 c
992 c Old arbitrary potential - commented out.
993 c
994 c      dbr= 4.20D0
995 c      fbr= 3.30D0
996 c
997 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
998 c energy surface of diethyl disulfide.
999 c A. Liwo and U. Kozlowska, 11/24/03
1000 c
1001       D0CM = 3.78d0
1002       AKCM = 15.1d0
1003       AKTH = 11.0d0
1004       AKCT = 12.0d0
1005       V1SS =-1.08d0
1006       V2SS = 7.61d0
1007       V3SS = 13.7d0
1008 c      akcm=0.0d0
1009 c      akth=0.0d0
1010 c      akct=0.0d0
1011 c      v1ss=0.0d0
1012 c      v2ss=0.0d0
1013 c      v3ss=0.0d0
1014       
1015       if(me.eq.king) then
1016       write (iout,'(/a)') "Disulfide bridge parameters:"
1017       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1018       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1019       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1020       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1021      &  ' v3ss:',v3ss
1022       endif
1023       return
1024   111 write (iout,*) "Error reading bending energy parameters."
1025       goto 999
1026   112 write (iout,*) "Error reading rotamer energy parameters."
1027       goto 999
1028   113 write (iout,*) "Error reading torsional energy parameters."
1029       goto 999
1030   114 write (iout,*) "Error reading double torsional energy parameters."
1031       goto 999
1032   115 write (iout,*) 
1033      &  "Error reading cumulant (multibody energy) parameters."
1034       goto 999
1035   116 write (iout,*) "Error reading electrostatic energy parameters."
1036       goto 999
1037   117 write (iout,*) "Error reading side chain interaction parameters."
1038       goto 999
1039   118 write (iout,*) "Error reading SCp interaction parameters."
1040       goto 999
1041   119 write (iout,*) "Error reading SCCOR parameters"
1042   999 continue
1043 #ifdef MPI
1044       call MPI_Finalize(Ierror)
1045 #endif
1046       stop
1047       return
1048       end
1049
1050
1051       subroutine getenv_loc(var, val)
1052       character(*) var, val
1053
1054 #ifdef WINIFL
1055       character(2000) line
1056       external ilen
1057
1058       open (196,file='env',status='old',readonly,shared)
1059       iread=0
1060 c      write(*,*)'looking for ',var
1061 10    read(196,*,err=11,end=11)line
1062       iread=index(line,var)
1063 c      write(*,*)iread,' ',var,' ',line
1064       if (iread.eq.0) go to 10 
1065 c      write(*,*)'---> ',line
1066 11    continue
1067       if(iread.eq.0) then
1068 c       write(*,*)'CHUJ'
1069        val=''
1070       else
1071        iread=iread+ilen(var)+1
1072        read (line(iread:),*,err=12,end=12) val
1073 c       write(*,*)'OK: ',var,' = ',val
1074       endif
1075       close(196)
1076       return
1077 12    val=''
1078       close(196)
1079 #elif (defined CRAY)
1080       integer lennam,lenval,ierror
1081 c
1082 c        getenv using a POSIX call, useful on the T3D
1083 c        Sept 1996, comment out error check on advice of H. Pritchard
1084 c
1085       lennam = len(var)
1086       if(lennam.le.0) stop '--error calling getenv--'
1087       call pxfgetenv(var,lennam,val,lenval,ierror)
1088 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1089 #else
1090       call getenv(var,val)
1091 #endif
1092
1093       return
1094       end