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