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