Merge branch 'master' of mmka.chem.univ.gda.pl:unres
[unres.git] / src-M / parmread.F
1       subroutine parmread(iparm,*)
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 'DIMENSIONS.ZSCOPT'
9       include 'DIMENSIONS.FREE'
10       include 'COMMON.IOUNITS'
11       include 'COMMON.CHAIN'
12       include 'COMMON.INTERACT'
13       include 'COMMON.GEO'
14       include 'COMMON.LOCAL'
15       include 'COMMON.TORSION'
16       include 'COMMON.FFIELD'
17       include 'COMMON.NAMES'
18       include 'COMMON.SBRIDGE'
19       include 'COMMON.WEIGHTS'
20       include 'COMMON.ENEPS'
21       include 'COMMON.SCCOR'
22       include 'COMMON.SCROT'
23       include 'COMMON.FREE'
24       character*1 t1,t2,t3
25       character*1 onelett(4) /"G","A","P","D"/
26       character*1 toronelet(-2:2) /"p","a","G","A","P"/
27       logical lprint
28       dimension blower(3,3,maxlob)
29       character*800 controlcard
30       character*256 bondname_t,thetname_t,rotname_t,torname_t,
31      & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
32      & sccorname_t
33       integer ilen
34       external ilen
35       character*16 key
36       integer iparm
37       double precision ip,mp
38       character*6 res1
39 C      write (iout,*) "KURWA"
40 C
41 C Body
42 C
43 C Set LPRINT=.TRUE. for debugging
44       dwa16=2.0d0**(1.0d0/6.0d0)
45       lprint=.false.
46       itypro=20
47 C Assign virtual-bond length
48       vbl=3.8D0
49       vblinv=1.0D0/vbl
50       vblinv2=vblinv*vblinv
51       call card_concat(controlcard,.true.)
52       wname(4)="WCORRH"
53       do i=1,n_ene
54         key = wname(i)(:ilen(wname(i)))
55         call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
56       enddo
57
58       write (iout,*) "iparm",iparm," myparm",myparm
59 c If reading not own parameters, skip assignment
60       call reada(controlcard,"D0CM",d0cm,3.78d0)
61       call reada(controlcard,"AKCM",akcm,15.1d0)
62       call reada(controlcard,"AKTH",akth,11.0d0)
63       call reada(controlcard,"AKCT",akct,12.0d0)
64       call reada(controlcard,"V1SS",v1ss,-1.08d0)
65       call reada(controlcard,"V2SS",v2ss,7.61d0)
66       call reada(controlcard,"V3SS",v3ss,13.7d0)
67       call reada(controlcard,"EBR",ebr,-5.50D0)
68       call reada(controlcard,"DTRISS",dtriss,1.0D0)
69       call reada(controlcard,"ATRISS",atriss,0.3D0)
70       call reada(controlcard,"BTRISS",btriss,0.02D0)
71       call reada(controlcard,"CTRISS",ctriss,1.0D0)
72       dyn_ss=(index(controlcard,'DYN_SS').gt.0)
73       write(iout,*) "ATRISS",atriss
74       write(iout,*) "BTRISS",btriss
75       write(iout,*) "CTRISS",ctriss
76       write(iout,*) "DTRISS",dtriss
77
78 C      do i=1,maxres
79 C        dyn_ss_mask(i)=.false.
80 C      enddo
81 C      ebr=-12.0D0
82 c
83 c Old arbitrary potential - commented out.
84 c
85 c      dbr= 4.20D0
86 c      fbr= 3.30D0
87 c
88 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
89 c energy surface of diethyl disulfide.
90 c A. Liwo and U. Kozlowska, 11/24/03
91 c
92       D0CM = 3.78d0
93       AKCM = 15.1d0
94       AKTH = 11.0d0
95       AKCT = 12.0d0
96       V1SS =-1.08d0
97       V2SS = 7.61d0
98       V3SS = 13.7d0
99
100       do i=1,maxres-1
101         do j=i+1,maxres
102           dyn_ssbond_ij(i,j)=1.0d300
103         enddo
104       enddo
105       call reada(controlcard,"HT",Ht,0.0D0)
106 C      if (dyn_ss) then
107 C        ss_depth=ebr/wsc-0.25*eps(1,1)
108 C        write(iout,*) HT,wsc,eps(1,1),'KURWA'
109 C        Ht=Ht/wsc-0.25*eps(1,1)
110        
111 C        akcm=akcm*whpb/wsc
112 C        akth=akth*whpb/wsc
113 C        akct=akct*whpb/wsc
114 C        v1ss=v1ss*whpb/wsc
115 C        v2ss=v2ss*whpb/wsc
116 C        v3ss=v3ss*whpb/wsc
117 C      else
118 C        ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
119 C      endif
120
121       if (iparm.eq.myparm .or. .not.separate_parset) then
122
123 c
124 c Setup weights for UNRES
125 c
126       wsc=ww(1)
127       wscp=ww(2)
128       welec=ww(3)
129       wcorr=ww(4)
130       wcorr5=ww(5)
131       wcorr6=ww(6)
132       wel_loc=ww(7)
133       wturn3=ww(8)
134       wturn4=ww(9)
135       wturn6=ww(10)
136       wang=ww(11)
137       wscloc=ww(12)
138       wtor=ww(13)
139       wtor_d=ww(14)
140       wvdwpp=ww(16)
141       wbond=ww(18)
142       wsccor=ww(19)
143       whpb=ww(15)
144       endif
145
146       call card_concat(controlcard,.false.)
147
148 c Return if not own parameters
149
150       if (iparm.ne.myparm .and. separate_parset) return
151
152       call reads(controlcard,"BONDPAR",bondname_t,bondname)
153       open (ibond,file=bondname_t,status='old')
154       rewind(ibond)
155       call reads(controlcard,"THETPAR",thetname_t,thetname)
156       open (ithep,file=thetname_t,status='old')
157       rewind(ithep) 
158       call reads(controlcard,"ROTPAR",rotname_t,rotname)
159       open (irotam,file=rotname_t,status='old')
160       rewind(irotam)
161       call reads(controlcard,"TORPAR",torname_t,torname)
162       open (itorp,file=torname_t,status='old')
163       rewind(itorp)
164       call reads(controlcard,"TORDPAR",tordname_t,tordname)
165       open (itordp,file=tordname_t,status='old')
166       rewind(itordp)
167       call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
168       open (isccor,file=sccorname_t,status='old')
169       rewind(isccor)
170       call reads(controlcard,"FOURIER",fouriername_t,fouriername)
171       open (ifourier,file=fouriername_t,status='old')
172       rewind(ifourier)
173       call reads(controlcard,"ELEPAR",elename_t,elename)
174       open (ielep,file=elename_t,status='old')
175       rewind(ielep)
176       call reads(controlcard,"SIDEPAR",sidename_t,sidename)
177       open (isidep,file=sidename_t,status='old')
178       rewind(isidep)
179       call reads(controlcard,"SCPPAR",scpname_t,scpname)
180       open (iscpp,file=scpname_t,status='old')
181       rewind(iscpp)
182       write (iout,*) "Parameter set:",iparm
183       write (iout,*) "Energy-term weights:"
184       do i=1,n_ene
185         write (iout,'(a16,f10.5)') wname(i),ww(i)
186       enddo
187       write (iout,*) "Sidechain potential file        : ",
188      &  sidename_t(:ilen(sidename_t))
189 #ifndef OLDSCP
190       write (iout,*) "SCp potential file              : ",
191      &  scpname_t(:ilen(scpname_t))
192 #endif  
193       write (iout,*) "Electrostatic potential file    : ",
194      &  elename_t(:ilen(elename_t))
195       write (iout,*) "Cumulant coefficient file       : ",
196      &  fouriername_t(:ilen(fouriername_t))
197       write (iout,*) "Torsional parameter file        : ",
198      &  torname_t(:ilen(torname_t))
199       write (iout,*) "Double torsional parameter file : ",
200      &  tordname_t(:ilen(tordname_t))
201       write (iout,*) "Backbone-rotamer parameter file : ",
202      &  sccorname(:ilen(sccorname))
203       write (iout,*) "Bond & inertia constant file    : ",
204      &  bondname_t(:ilen(bondname_t))
205       write (iout,*) "Bending parameter file          : ",
206      &  thetname_t(:ilen(thetname_t))
207       write (iout,*) "Rotamer parameter file          : ",
208      &  rotname_t(:ilen(rotname_t))
209
210 c
211 c Read the virtual-bond parameters, masses, and moments of inertia
212 c and Stokes' radii of the peptide group and side chains
213 c
214 #ifdef CRYST_BOND
215       read (ibond,*) vbldp0,akp
216       do i=1,ntyp
217         nbondterm(i)=1
218         read (ibond,*) vbldsc0(1,i),aksc(1,i)
219         dsc(i) = vbldsc0(1,i)
220         if (i.eq.10) then
221           dsc_inv(i)=0.0D0
222         else
223           dsc_inv(i)=1.0D0/dsc(i)
224         endif
225       enddo
226 #else
227       read (ibond,*) ijunk,vbldp0,akp,rjunk
228       do i=1,ntyp
229         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
230      &   j=1,nbondterm(i))
231         dsc(i) = vbldsc0(1,i)
232         if (i.eq.10) then
233           dsc_inv(i)=0.0D0
234         else
235           dsc_inv(i)=1.0D0/dsc(i)
236         endif
237       enddo
238 #endif
239       if (lprint) then
240         write(iout,'(/a/)')"Force constants virtual bonds:"
241         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
242      &   'inertia','Pstok'
243         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
244         do i=1,ntyp
245           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
246      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
247           do j=2,nbondterm(i)
248             write (iout,'(13x,3f10.5)')
249      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
250           enddo
251         enddo
252       endif
253 #ifdef CRYST_THETA
254 C
255 C Read the parameters of the probability distribution/energy expression 
256 C of the virtual-bond valence angles theta
257 C
258       do i=1,ntyp
259         read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
260      &    (bthet(j,i,1,1),j=1,2)
261         read (ithep,*) (polthet(j,i),j=0,3)
262         read (ithep,*) (gthet(j,i),j=1,3)
263         read (ithep,*) theta0(i),sig0(i),sigc0(i)
264         sigc0(i)=sigc0(i)**2
265       enddo
266       do i=1,ntyp
267       athet(1,i,1,-1)=athet(1,i,1,1)
268       athet(2,i,1,-1)=athet(2,i,1,1)
269       bthet(1,i,1,-1)=-bthet(1,i,1,1)
270       bthet(2,i,1,-1)=-bthet(2,i,1,1)
271       athet(1,i,-1,1)=-athet(1,i,1,1)
272       athet(2,i,-1,1)=-athet(2,i,1,1)
273       bthet(1,i,-1,1)=bthet(1,i,1,1)
274       bthet(2,i,-1,1)=bthet(2,i,1,1)
275       enddo
276       do i=-ntyp,-1
277       a0thet(i)=a0thet(-i)
278       athet(1,i,-1,-1)=athet(1,-i,1,1)
279       athet(2,i,-1,-1)=-athet(2,-i,1,1)
280       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
281       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
282       athet(1,i,-1,1)=athet(1,-i,1,1)
283       athet(2,i,-1,1)=-athet(2,-i,1,1)
284       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
285       bthet(2,i,-1,1)=bthet(2,-i,1,1)
286       athet(1,i,1,-1)=-athet(1,-i,1,1)
287       athet(2,i,1,-1)=athet(2,-i,1,1)
288       bthet(1,i,1,-1)=bthet(1,-i,1,1)
289       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
290       theta0(i)=theta0(-i)
291       sig0(i)=sig0(-i)
292       sigc0(i)=sigc0(-i)
293        do j=0,3
294         polthet(j,i)=polthet(j,-i)
295        enddo
296        do j=1,3
297          gthet(j,i)=gthet(j,-i)
298        enddo
299       enddo
300       close (ithep)
301       if (lprint) then
302 c       write (iout,'(a)') 
303 c    &    'Parameters of the virtual-bond valence angles:'
304 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
305 c    & '    ATHETA0   ','         A1   ','        A2    ',
306 c    & '        B1    ','         B2   '        
307 c       do i=1,ntyp
308 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
309 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
310 c       enddo
311 c       write (iout,'(/a/9x,5a/79(1h-))') 
312 c    & 'Parameters of the expression for sigma(theta_c):',
313 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
314 c    & '     ALPH3    ','    SIGMA0C   '        
315 c       do i=1,ntyp
316 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
317 c    &      (polthet(j,i),j=0,3),sigc0(i) 
318 c       enddo
319 c       write (iout,'(/a/9x,5a/79(1h-))') 
320 c    & 'Parameters of the second gaussian:',
321 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
322 c    & '        G2    ','         G3   '        
323 c       do i=1,ntyp
324 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
325 c    &       sig0(i),(gthet(j,i),j=1,3)
326 c       enddo
327         write (iout,'(a)') 
328      &    'Parameters of the virtual-bond valence angles:'
329         write (iout,'(/a/9x,5a/79(1h-))') 
330      & 'Coefficients of expansion',
331      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
332      & '   b1*10^1    ','    b2*10^1   '        
333         do i=1,ntyp
334           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
335      &        a0thet(i),(100*athet(j,i,1,1),j=1,2),
336      &        (10*bthet(j,i,1,1),j=1,2)
337         enddo
338         write (iout,'(/a/9x,5a/79(1h-))') 
339      & 'Parameters of the expression for sigma(theta_c):',
340      & ' alpha0       ','  alph1       ',' alph2        ',
341      & ' alhp3        ','   sigma0c    '        
342         do i=1,ntyp
343           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
344      &      (polthet(j,i),j=0,3),sigc0(i) 
345         enddo
346         write (iout,'(/a/9x,5a/79(1h-))') 
347      & 'Parameters of the second gaussian:',
348      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
349      & '        G2    ','   G3*10^1    '        
350         do i=1,ntyp
351           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
352      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
353         enddo
354       endif
355 #else
356 C
357 C Read the parameters of Utheta determined from ab initio surfaces
358 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
359 C
360       write (iout,*) "tu dochodze"
361       read (ithep,*) nthetyp,ntheterm,ntheterm2,
362      &  ntheterm3,nsingle,ndouble
363       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
364       read (ithep,*) (ithetyp(i),i=1,ntyp1)
365       do i=-ntyp1,-1
366         ithetyp(i)=-ithetyp(-i)
367       enddo
368       write (iout,*) "tu dochodze"
369       do iblock=1,2
370       do i=-maxthetyp,maxthetyp
371         do j=-maxthetyp,maxthetyp
372           do k=-maxthetyp,maxthetyp
373             aa0thet(i,j,k,iblock)=0.0d0
374             do l=1,ntheterm
375               aathet(l,i,j,k,iblock)=0.0d0
376             enddo
377             do l=1,ntheterm2
378               do m=1,nsingle
379                 bbthet(m,l,i,j,k,iblock)=0.0d0
380                 ccthet(m,l,i,j,k,iblock)=0.0d0
381                 ddthet(m,l,i,j,k,iblock)=0.0d0
382                 eethet(m,l,i,j,k,iblock)=0.0d0
383               enddo
384             enddo
385             do l=1,ntheterm3
386               do m=1,ndouble
387                 do mm=1,ndouble
388                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
389                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
390                 enddo
391               enddo
392             enddo
393           enddo
394         enddo
395       enddo
396       enddo
397 C      write (iout,*) "KURWA1"
398       do iblock=1,2
399       do i=0,nthetyp
400         do j=-nthetyp,nthetyp
401           do k=-nthetyp,nthetyp
402             read (ithep,'(6a)') res1
403             write(iout,*) res1,i,j,k
404             read (ithep,*) aa0thet(i,j,k,iblock)
405             read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
406             read (ithep,*)
407      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
408      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
409      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
410      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
411      &        ,ll=1,ntheterm2)
412             read (ithep,*)
413      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
414      &      ffthet(lll,llll,ll,i,j,k,iblock),
415      &         ggthet(llll,lll,ll,i,j,k,iblock)
416      &        ,ggthet(lll,llll,ll,i,j,k,iblock),
417      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
418           enddo
419         enddo
420       enddo
421 C       write(iout,*) "KURWA1.1"
422 C
423 C For dummy ends assign glycine-type coefficients of theta-only terms; the
424 C coefficients of theta-and-gamma-dependent terms are zero.
425 C
426       do i=1,nthetyp
427         do j=1,nthetyp
428           do l=1,ntheterm
429             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
430             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
431           enddo
432           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
433           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
434         enddo
435         do l=1,ntheterm
436           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
437         enddo
438         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
439       enddo
440       enddo
441 C       write(iout,*) "KURWA1.5"
442 C Substitution for D aminoacids from symmetry.
443       do iblock=1,2
444       do i=-nthetyp,0
445         do j=-nthetyp,nthetyp
446           do k=-nthetyp,nthetyp
447            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
448            do l=1,ntheterm
449            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
450            enddo
451            do ll=1,ntheterm2
452             do lll=1,nsingle
453             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
454             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
455             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
456             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
457             enddo
458           enddo
459           do ll=1,ntheterm3
460            do lll=2,ndouble
461             do llll=1,lll-1
462             ffthet(llll,lll,ll,i,j,k,iblock)=
463      &      ffthet(llll,lll,ll,-i,-j,-k,iblock)
464             ffthet(lll,llll,ll,i,j,k,iblock)=
465      &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
466             ggthet(llll,lll,ll,i,j,k,iblock)=
467      &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
468             ggthet(lll,llll,ll,i,j,k,iblock)=
469      &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)
470             enddo !ll
471            enddo  !lll  
472           enddo   !llll
473          enddo    !k
474         enddo     !j
475        enddo      !i
476       enddo       !iblock
477
478 C
479 C Control printout of the coefficients of virtual-bond-angle potentials
480 C
481       if (lprint) then
482         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
483         do i=1,nthetyp+1
484           do j=1,nthetyp+1
485             do k=1,nthetyp+1
486               write (iout,'(//4a)')
487      &         'Type ',onelett(i),onelett(j),onelett(k)
488               write (iout,'(//a,10x,a)') " l","a[l]"
489               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
490               write (iout,'(i2,1pe15.5)')
491      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
492             do l=1,ntheterm2
493               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
494      &          "b",l,"c",l,"d",l,"e",l
495               do m=1,nsingle
496                 write (iout,'(i2,4(1pe15.5))') m,
497      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
498      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
499               enddo
500             enddo
501             do l=1,ntheterm3
502               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
503      &          "f+",l,"f-",l,"g+",l,"g-",l
504               do m=2,ndouble
505                 do n=1,m-1
506                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
507      &              ffthet(n,m,l,i,j,k,iblock),
508      &              ffthet(m,n,l,i,j,k,iblock),
509      &              ggthet(n,m,l,i,j,k,iblock),
510      &              ggthet(m,n,l,i,j,k,iblock)
511                 enddo
512               enddo
513             enddo
514           enddo
515         enddo
516       enddo
517       call flush(iout)
518       endif
519 #endif
520 C      write(iout,*) 'KURWA2'
521 #ifdef CRYST_SC
522 C
523 C Read the parameters of the probability distribution/energy expression
524 C of the side chains.
525 C
526       do i=1,ntyp
527 cc      write (iout,*) "tu dochodze",i
528         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
529         if (i.eq.10) then
530           dsc_inv(i)=0.0D0
531         else
532           dsc_inv(i)=1.0D0/dsc(i)
533         endif
534         if (i.ne.10) then
535         do j=1,nlob(i)
536           do k=1,3
537             do l=1,3
538               blower(l,k,j)=0.0D0
539             enddo
540           enddo
541         enddo  
542         bsc(1,i)=0.0D0
543         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
544         censc(1,1,-i)=censc(1,1,i)
545         censc(2,1,-i)=censc(2,1,i)
546         censc(3,1,-i)=-censc(3,1,i)
547         do j=2,nlob(i)
548           read (irotam,*) bsc(j,i)
549           read (irotam,*) (censc(k,j,i),k=1,3),
550      &                                 ((blower(k,l,j),l=1,k),k=1,3)
551         censc(1,j,-i)=censc(1,j,i)
552         censc(2,j,-i)=censc(2,j,i)
553         censc(3,j,-i)=-censc(3,j,i)
554 C BSC is amplitude of Gaussian
555         enddo
556         do j=1,nlob(i)
557           do k=1,3
558             do l=1,k
559               akl=0.0D0
560               do m=1,3
561                 akl=akl+blower(k,m,j)*blower(l,m,j)
562               enddo
563               gaussc(k,l,j,i)=akl
564               gaussc(l,k,j,i)=akl
565              if (((k.eq.3).and.(l.ne.3))
566      &        .or.((l.eq.3).and.(k.ne.3))) then
567                 gaussc(k,l,j,-i)=-akl
568                 gaussc(l,k,j,-i)=-akl
569               else
570                 gaussc(k,l,j,-i)=akl
571                 gaussc(l,k,j,-i)=akl
572               endif
573             enddo
574           enddo 
575         enddo
576         endif
577       enddo
578       close (irotam)
579       if (lprint) then
580         write (iout,'(/a)') 'Parameters of side-chain local geometry'
581         do i=1,ntyp
582           nlobi=nlob(i)
583           if (nlobi.gt.0) then
584           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
585      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
586 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
587 c          write (iout,'(a,f10.4,4(16x,f10.4))')
588 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
589 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
590            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
591      &                             'log h',(bsc(j,i),j=1,nlobi)
592            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
593      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
594 c          write (iout,'(a)')
595 c         do j=1,nlobi
596 c           ind=0
597 c           do k=1,3
598 c             do l=1,k
599 c              ind=ind+1
600 c              blower(k,l,j)=gaussc(ind,j,i)
601 c             enddo
602 c           enddo
603 c         enddo
604           do k=1,3
605             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
606      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
607           enddo
608           endif
609         enddo
610       endif
611 #else
612 C
613 C Read scrot parameters for potentials determined from all-atom AM1 calculations
614 C added by Urszula Kozlowska 07/11/2007
615 C
616       do i=1,ntyp
617         read (irotam,*)
618        if (i.eq.10) then
619          read (irotam,*)
620        else
621          do j=1,65
622            read(irotam,*) sc_parmin(j,i)
623          enddo
624        endif
625       enddo
626 #endif
627       close(irotam)
628 C      write (iout,*) 'KURWAKURWA'
629 #ifdef CRYST_TOR
630 C
631 C Read torsional parameters in old format
632 C
633       read (itorp,*) ntortyp,nterm_old
634       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
635       read (itorp,*) (itortyp(i),i=1,ntyp)
636       do i=1,ntortyp
637         do j=1,ntortyp
638           read (itorp,'(a)')
639           do k=1,nterm_old
640             read (itorp,*) kk,v1(k,j,i),v2(k,j,i) 
641           enddo
642         enddo
643       enddo
644       close (itorp)
645       if (lprint) then
646         write (iout,'(/a/)') 'Torsional constants:'
647         do i=1,ntortyp
648           do j=1,ntortyp
649             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
650             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
651           enddo
652         enddo
653       endif
654
655
656 #else
657 C
658 C Read torsional parameters
659 C
660       read (itorp,*) ntortyp
661       read (itorp,*) (itortyp(i),i=1,ntyp)
662       write (iout,*) 'ntortyp',ntortyp
663       do iblock=1,2
664       do i=-ntyp,-1
665        itortyp(i)=-itortyp(-i)
666       enddo
667 c      write (iout,*) 'ntortyp',ntortyp
668       do i=0,ntortyp-1
669         do j=-ntortyp+1,ntortyp-1
670           read (itorp,*) nterm(i,j,iblock),
671      &          nlor(i,j,iblock)
672           nterm(-i,-j,iblock)=nterm(i,j,iblock)
673           nlor(-i,-j,iblock)=nlor(i,j,iblock)
674           v0ij=0.0d0
675           si=-1.0d0
676           do k=1,nterm(i,j,iblock)
677             read (itorp,*) kk,v1(k,i,j,iblock),
678      &      v2(k,i,j,iblock)
679             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
680             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
681             v0ij=v0ij+si*v1(k,i,j,iblock)
682             si=-si
683          enddo
684           do k=1,nlor(i,j,iblock)
685             read (itorp,*) kk,vlor1(k,i,j),
686      &        vlor2(k,i,j),vlor3(k,i,j)
687             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
688           enddo
689           v0(i,j,iblock)=v0ij
690           v0(-i,-j,iblock)=v0ij
691         enddo
692       enddo
693       enddo
694       close (itorp)
695       if (lprint) then
696         write (iout,'(/a/)') 'Torsional constants:'
697         do i=1,ntortyp
698           do j=1,ntortyp
699             write (iout,*) 'ityp',i,' jtyp',j
700             write (iout,*) 'Fourier constants'
701             do k=1,nterm(i,j,iblock)
702               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
703      &        v2(k,i,j,iblock)
704             enddo
705             write (iout,*) 'Lorenz constants'
706             do k=1,nlor(i,j,iblock)
707               write (iout,'(3(1pe15.5))')
708      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
709             enddo
710           enddo
711         enddo
712       endif
713 C
714 C 6/23/01 Read parameters for double torsionals
715 C
716       do iblock=1,2
717       do i=0,ntortyp-1
718         do j=-ntortyp+1,ntortyp-1
719           do k=-ntortyp+1,ntortyp-1
720             read (itordp,'(3a1)') t1,t2,t3
721 c              write (iout,*) "OK onelett",
722 c     &         i,j,k,t1,t2,t3
723
724             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
725      &        .or. t3.ne.toronelet(k)) then
726               write (iout,*) "Error in double torsional parameter file",
727      &         i,j,k,t1,t2,t3
728 #ifdef MPI
729               call MPI_Finalize(Ierror)
730 #endif
731                stop "Error in double torsional parameter file"
732             endif
733           read (itordp,*) ntermd_1(i,j,k,iblock),
734      &         ntermd_2(i,j,k,iblock)
735             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
736             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
737             read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
738      &         ntermd_1(i,j,k,iblock))
739             read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
740      &         ntermd_1(i,j,k,iblock))
741             read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
742      &         ntermd_1(i,j,k,iblock))
743             read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
744      &         ntermd_1(i,j,k,iblock))
745 C Martix of D parameters for one dimesional foureir series
746             do l=1,ntermd_1(i,j,k,iblock)
747              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
748              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
749              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
750              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
751 c            write(iout,*) "whcodze" ,
752 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
753             enddo
754             read (itordp,*) ((v2c(l,m,i,j,k,iblock),
755      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
756      &         v2s(m,l,i,j,k,iblock),
757      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
758 C Martix of D parameters for two dimesional fourier series
759             do l=1,ntermd_2(i,j,k,iblock)
760              do m=1,l-1
761              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
762              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
763              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
764              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
765              enddo!m
766             enddo!l
767           enddo!k
768         enddo!j
769       enddo!i
770       enddo!iblock
771       if (lprint) then
772       write (iout,*)
773       write (iout,*) 'Constants for double torsionals'
774       do iblock=1,2
775       do i=0,ntortyp-1
776         do j=-ntortyp+1,ntortyp-1
777           do k=-ntortyp+1,ntortyp-1
778             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
779      &        ' nsingle',ntermd_1(i,j,k,iblock),
780      &        ' ndouble',ntermd_2(i,j,k,iblock)
781             write (iout,*)
782             write (iout,*) 'Single angles:'
783             do l=1,ntermd_1(i,j,k,iblock)
784               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
785      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
786      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
787      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
788             enddo
789             write (iout,*)
790             write (iout,*) 'Pairs of angles:'
791             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
792             do l=1,ntermd_2(i,j,k,iblock)
793               write (iout,'(i5,20f10.5)')
794      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
795             enddo
796             write (iout,*)
797            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
798             do l=1,ntermd_2(i,j,k,iblock)
799               write (iout,'(i5,20f10.5)')
800      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
801      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
802             enddo
803             write (iout,*)
804           enddo
805         enddo
806       enddo
807       enddo
808       endif
809 #endif
810 C Read of Side-chain backbone correlation parameters
811 C Modified 11 May 2012 by Adasko
812 CCC
813 C
814       read (isccor,*) nsccortyp
815       read (isccor,*) (isccortyp(i),i=1,ntyp)
816       do i=-ntyp,-1
817         isccortyp(i)=-isccortyp(-i)
818       enddo
819       iscprol=isccortyp(20)
820 c      write (iout,*) 'ntortyp',ntortyp
821       maxinter=3
822 cc maxinter is maximum interaction sites
823       do l=1,maxinter
824       do i=1,nsccortyp
825         do j=1,nsccortyp
826           read (isccor,*)
827      &nterm_sccor(i,j),nlor_sccor(i,j)
828           write (iout,*) nterm_sccor(i,j)
829           v0ijsccor=0.0d0
830           v0ijsccor1=0.0d0
831           v0ijsccor2=0.0d0
832           v0ijsccor3=0.0d0
833           si=-1.0d0
834           nterm_sccor(-i,j)=nterm_sccor(i,j)
835           nterm_sccor(-i,-j)=nterm_sccor(i,j)
836           nterm_sccor(i,-j)=nterm_sccor(i,j)
837           write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
838      &    nterm_sccor(-i,-j),nterm_sccor(i,-j)
839           do k=1,nterm_sccor(i,j)
840             read (isccor,*) kk,v1sccor(k,l,i,j)
841      &    ,v2sccor(k,l,i,j)
842             if (j.eq.iscprol) then
843              if (i.eq.isccortyp(10)) then
844              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
845              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
846              else
847              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
848      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
849              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
850      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
851              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
852              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
853              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
854              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
855              endif
856             else
857              if (i.eq.isccortyp(10)) then
858              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
859              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
860              else
861                if (j.eq.isccortyp(10)) then
862              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
863              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
864                else
865              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
866              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
867              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
868              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
869              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
870              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
871                 endif
872                endif
873             endif
874             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
875             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
876             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
877             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
878             si=-si
879            enddo
880           do k=1,nlor_sccor(i,j)
881             read (isccor,*) kk,vlor1sccor(k,i,j),
882      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
883             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
884      &(1+vlor3sccor(k,i,j)**2)
885           enddo
886           v0sccor(l,i,j)=v0ijsccor
887           v0sccor(l,-i,j)=v0ijsccor1
888           v0sccor(l,i,-j)=v0ijsccor2
889           v0sccor(l,-i,-j)=v0ijsccor3
890           enddo
891         enddo
892       enddo
893       close (isccor)
894       if (lprint) then
895         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
896         do i=1,nsccortyp
897           do j=1,nsccortyp
898             write (iout,*) 'ityp',i,' jtyp',j
899             write (iout,*) 'Fourier constants'
900             do k=1,nterm_sccor(i,j)
901               write (iout,'(2(1pe15.5))')
902      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
903             enddo
904             write (iout,*) 'Lorenz constants'
905             do k=1,nlor_sccor(i,j)
906               write (iout,'(3(1pe15.5))')
907      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
908             enddo
909           enddo
910         enddo
911       endif
912 C
913 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
914 C         interaction energy of the Gly, Ala, and Pro prototypes.
915 C
916       read (ifourier,*) nloctyp
917       do i=0,nloctyp-1
918         read (ifourier,*)
919         read (ifourier,*) (b(ii,i),ii=1,13)
920         if (lprint) then
921         write (iout,*) 'Type',i
922         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
923         endif
924         B1(1,i)  = b(3,i)
925         B1(2,i)  = b(5,i)
926         B1(1,-i) = b(3,i)
927         B1(2,-i) = -b(5,i)
928 c        b1(1,i)=0.0d0
929 c        b1(2,i)=0.0d0
930         B1tilde(1,i) = b(3,i)
931         B1tilde(2,i) =-b(5,i)
932         B1tilde(1,-i) =-b(3,i)
933         B1tilde(2,-i) =b(5,i)
934 c        b1tilde(1,i)=0.0d0
935 c        b1tilde(2,i)=0.0d0
936         B2(1,i)  = b(2,i)
937         B2(2,i)  = b(4,i)
938         B2(1,-i)  =b(2,i)
939         B2(2,-i)  =-b(4,i)
940
941 c        b2(1,i)=0.0d0
942 c        b2(2,i)=0.0d0
943         CC(1,1,i)= b(7,i)
944         CC(2,2,i)=-b(7,i)
945         CC(2,1,i)= b(9,i)
946         CC(1,2,i)= b(9,i)
947         CC(1,1,-i)= b(7,i)
948         CC(2,2,-i)=-b(7,i)
949         CC(2,1,-i)=-b(9,i)
950         CC(1,2,-i)=-b(9,i)
951 c        CC(1,1,i)=0.0d0
952 c        CC(2,2,i)=0.0d0
953 c        CC(2,1,i)=0.0d0
954 c        CC(1,2,i)=0.0d0
955         Ctilde(1,1,i)=b(7,i)
956         Ctilde(1,2,i)=b(9,i)
957         Ctilde(2,1,i)=-b(9,i)
958         Ctilde(2,2,i)=b(7,i)
959         Ctilde(1,1,-i)=b(7,i)
960         Ctilde(1,2,-i)=-b(9,i)
961         Ctilde(2,1,-i)=b(9,i)
962         Ctilde(2,2,-i)=b(7,i)
963
964 c        Ctilde(1,1,i)=0.0d0
965 c        Ctilde(1,2,i)=0.0d0
966 c        Ctilde(2,1,i)=0.0d0
967 c        Ctilde(2,2,i)=0.0d0
968         DD(1,1,i)= b(6,i)
969         DD(2,2,i)=-b(6,i)
970         DD(2,1,i)= b(8,i)
971         DD(1,2,i)= b(8,i)
972         DD(1,1,-i)= b(6,i)
973         DD(2,2,-i)=-b(6,i)
974         DD(2,1,-i)=-b(8,i)
975         DD(1,2,-i)=-b(8,i)
976 c        DD(1,1,i)=0.0d0
977 c        DD(2,2,i)=0.0d0
978 c        DD(2,1,i)=0.0d0
979 c        DD(1,2,i)=0.0d0
980         Dtilde(1,1,i)=b(6,i)
981         Dtilde(1,2,i)=b(8,i)
982         Dtilde(2,1,i)=-b(8,i)
983         Dtilde(2,2,i)=b(6,i)
984         Dtilde(1,1,-i)=b(6,i)
985         Dtilde(1,2,-i)=-b(8,i)
986         Dtilde(2,1,-i)=b(8,i)
987         Dtilde(2,2,-i)=b(6,i)
988
989 c        Dtilde(1,1,i)=0.0d0
990 c        Dtilde(1,2,i)=0.0d0
991 c        Dtilde(2,1,i)=0.0d0
992 c        Dtilde(2,2,i)=0.0d0
993         EE(1,1,i)= b(10,i)+b(11,i)
994         EE(2,2,i)=-b(10,i)+b(11,i)
995         EE(2,1,i)= b(12,i)-b(13,i)
996         EE(1,2,i)= b(12,i)+b(13,i)
997         EE(1,1,-i)= b(10,i)+b(11,i)
998         EE(2,2,-i)=-b(10,i)+b(11,i)
999         EE(2,1,-i)=-b(12,i)+b(13,i)
1000         EE(1,2,-i)=-b(12,i)-b(13,i)
1001
1002 c        ee(1,1,i)=1.0d0
1003 c        ee(2,2,i)=1.0d0
1004 c        ee(2,1,i)=0.0d0
1005 c        ee(1,2,i)=0.0d0
1006 c        ee(2,1,i)=ee(1,2,i)
1007
1008       enddo
1009       if (lprint) then
1010       do i=1,nloctyp
1011         write (iout,*) 'Type',i
1012         write (iout,*) 'B1'
1013 c        write (iout,'(f10.5)') B1(:,i)
1014         write(iout,*) B1(1,i),B1(2,i)
1015         write (iout,*) 'B2'
1016 c        write (iout,'(f10.5)') B2(:,i)
1017         write(iout,*) B2(1,i),B2(2,i)
1018         write (iout,*) 'CC'
1019         do j=1,2
1020           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1021         enddo
1022         write(iout,*) 'DD'
1023         do j=1,2
1024           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1025         enddo
1026         write(iout,*) 'EE'
1027         do j=1,2
1028           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1029         enddo
1030       enddo
1031       endif
1032
1033 C Read electrostatic-interaction parameters
1034 C
1035       if (lprint) then
1036         write (iout,'(/a)') 'Electrostatic interaction constants:'
1037         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
1038      &            'IT','JT','APP','BPP','AEL6','AEL3'
1039       endif
1040       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1041       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1042       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1043       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1044       close (ielep)
1045       do i=1,2
1046         do j=1,2
1047         rri=rpp(i,j)**6
1048         app (i,j)=epp(i,j)*rri*rri 
1049         bpp (i,j)=-2.0D0*epp(i,j)*rri
1050         ael6(i,j)=elpp6(i,j)*4.2D0**6
1051         ael3(i,j)=elpp3(i,j)*4.2D0**3
1052         lprint=.true.
1053         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1054      &                    ael6(i,j),ael3(i,j)
1055          lprint=.false.
1056         enddo
1057       enddo
1058 C
1059 C Read side-chain interaction parameters.
1060 C
1061       read (isidep,*) ipot,expon
1062       if (ipot.lt.1 .or. ipot.gt.5) then
1063         write (iout,'(2a)') 'Error while reading SC interaction',
1064      &               'potential file - unknown potential type.'
1065         stop
1066       endif
1067       expon2=expon/2
1068       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1069      & ', exponents are ',expon,2*expon 
1070       goto (10,20,30,30,40) ipot
1071 C----------------------- LJ potential ---------------------------------
1072    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1073       if (lprint) then
1074         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1075         write (iout,'(a/)') 'The epsilon array:'
1076         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1077         write (iout,'(/a)') 'One-body parameters:'
1078         write (iout,'(a,4x,a)') 'residue','sigma'
1079         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1080       endif
1081       goto 50
1082 C----------------------- LJK potential --------------------------------
1083    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1084      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1085       if (lprint) then
1086         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1087         write (iout,'(a/)') 'The epsilon array:'
1088         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1089         write (iout,'(/a)') 'One-body parameters:'
1090         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1091         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1092      &        i=1,ntyp)
1093       endif
1094       goto 50
1095 C---------------------- GB or BP potential -----------------------------
1096    30 do i=1,ntyp
1097        read (isidep,*)(eps(i,j),j=i,ntyp)
1098       enddo
1099       read (isidep,*)(sigma0(i),i=1,ntyp)
1100       read (isidep,*)(sigii(i),i=1,ntyp)
1101       read (isidep,*)(chip(i),i=1,ntyp)
1102       read (isidep,*)(alp(i),i=1,ntyp)
1103 C For the GB potential convert sigma'**2 into chi'
1104       if (ipot.eq.4) then
1105         do i=1,ntyp
1106           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1107         enddo
1108       endif
1109       if (lprint) then
1110         write (iout,'(/a/)') 'Parameters of the BP potential:'
1111         write (iout,'(a/)') 'The epsilon array:'
1112         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1113         write (iout,'(/a)') 'One-body parameters:'
1114         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1115      &       '    chip  ','    alph  '
1116         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1117      &                     chip(i),alp(i),i=1,ntyp)
1118       endif
1119       goto 50
1120 C--------------------- GBV potential -----------------------------------
1121    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1122      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1123      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1124       if (lprint) then
1125         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1126         write (iout,'(a/)') 'The epsilon array:'
1127         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1128         write (iout,'(/a)') 'One-body parameters:'
1129         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1130      &      's||/s_|_^2','    chip  ','    alph  '
1131         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1132      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1133       endif
1134    50 continue
1135       close (isidep)
1136 C-----------------------------------------------------------------------
1137 C Calculate the "working" parameters of SC interactions.
1138       do i=2,ntyp
1139         do j=1,i-1
1140           eps(i,j)=eps(j,i)
1141         enddo
1142       enddo
1143       do i=1,ntyp
1144         do j=i,ntyp
1145           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1146           sigma(j,i)=sigma(i,j)
1147           rs0(i,j)=dwa16*sigma(i,j)
1148           rs0(j,i)=rs0(i,j)
1149         enddo
1150       enddo
1151       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
1152      & 'Working parameters of the SC interactions:',
1153      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1154      & '  chi1   ','   chi2   ' 
1155       do i=1,ntyp
1156         do j=i,ntyp
1157           epsij=eps(i,j)
1158           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1159             rrij=sigma(i,j)
1160           else
1161             rrij=rr0(i)+rr0(j)
1162           endif
1163           r0(i,j)=rrij
1164           r0(j,i)=rrij
1165           rrij=rrij**expon
1166           epsij=eps(i,j)
1167           sigeps=dsign(1.0D0,epsij)
1168           epsij=dabs(epsij)
1169           aa(i,j)=epsij*rrij*rrij
1170           bb(i,j)=-sigeps*epsij*rrij
1171           aa(j,i)=aa(i,j)
1172           bb(j,i)=bb(i,j)
1173           if (ipot.gt.2) then
1174             sigt1sq=sigma0(i)**2
1175             sigt2sq=sigma0(j)**2
1176             sigii1=sigii(i)
1177             sigii2=sigii(j)
1178             ratsig1=sigt2sq/sigt1sq
1179             ratsig2=1.0D0/ratsig1
1180             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1181             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1182             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1183           else
1184             rsum_max=sigma(i,j)
1185           endif
1186 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1187             sigmaii(i,j)=rsum_max
1188             sigmaii(j,i)=rsum_max 
1189 c         else
1190 c           sigmaii(i,j)=r0(i,j)
1191 c           sigmaii(j,i)=r0(i,j)
1192 c         endif
1193 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1194           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1195             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1196             augm(i,j)=epsij*r_augm**(2*expon)
1197 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1198             augm(j,i)=augm(i,j)
1199           else
1200             augm(i,j)=0.0D0
1201             augm(j,i)=0.0D0
1202           endif
1203           if (lprint) then
1204             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
1205      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1206      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1207           endif
1208         enddo
1209       enddo
1210 C
1211 C Define the SC-p interaction constants
1212 C
1213 #ifdef OLDSCP
1214       do i=1,20
1215 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates 
1216 C helix formation)
1217 c       aad(i,1)=0.3D0*4.0D0**12
1218 C Following line for constants currently implemented
1219 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1220         aad(i,1)=1.5D0*4.0D0**12
1221 c       aad(i,1)=0.17D0*5.6D0**12
1222         aad(i,2)=aad(i,1)
1223 C "Soft" SC-p repulsion
1224         bad(i,1)=0.0D0
1225 C Following line for constants currently implemented
1226 c       aad(i,1)=0.3D0*4.0D0**6
1227 C "Hard" SC-p repulsion
1228         bad(i,1)=3.0D0*4.0D0**6
1229 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1230         bad(i,2)=bad(i,1)
1231 c       aad(i,1)=0.0D0
1232 c       aad(i,2)=0.0D0
1233 c       bad(i,1)=1228.8D0
1234 c       bad(i,2)=1228.8D0
1235       enddo
1236 #else
1237 C
1238 C 8/9/01 Read the SC-p interaction constants from file
1239 C
1240       do i=1,ntyp
1241         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1242       enddo
1243       do i=1,ntyp
1244         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1245         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1246         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1247         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1248       enddo
1249
1250       if (lprint) then
1251         write (iout,*) "Parameters of SC-p interactions:"
1252         do i=1,20
1253           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1254      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1255         enddo
1256       endif
1257 #endif
1258 C
1259 C Define the constants of the disulfide bridge
1260 C
1261 C      ebr=-12.0D0
1262 c
1263 c Old arbitrary potential - commented out.
1264 c
1265 c      dbr= 4.20D0
1266 c      fbr= 3.30D0
1267 c
1268 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1269 c energy surface of diethyl disulfide.
1270 c A. Liwo and U. Kozlowska, 11/24/03
1271 c
1272 C      D0CM = 3.78d0
1273 C      AKCM = 15.1d0
1274 C      AKTH = 11.0d0
1275 C      AKCT = 12.0d0
1276 C      V1SS =-1.08d0
1277 C      V2SS = 7.61d0
1278 C      V3SS = 13.7d0
1279       write (iout,*) dyn_ss,'dyndyn'
1280       if (dyn_ss) then
1281         ss_depth=ebr/wsc-0.25*eps(1,1)
1282 C        write(iout,*) akcm,whpb,wsc,'KURWA'
1283         Ht=Ht/wsc-0.25*eps(1,1)
1284
1285         akcm=akcm*whpb/wsc
1286         akth=akth*whpb/wsc
1287         akct=akct*whpb/wsc
1288         v1ss=v1ss*whpb/wsc
1289         v2ss=v2ss*whpb/wsc
1290         v3ss=v3ss*whpb/wsc
1291       else
1292         ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1293       endif
1294
1295 C      if (lprint) then
1296       write (iout,'(/a)') "Disulfide bridge parameters:"
1297       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1298       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1299       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1300       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1301      & ' v3ss:',v3ss
1302 C      endif
1303       return
1304       end