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