8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX,SPLIT_FOURIERTOR
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj
743 integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
749 integer :: ichir1,ichir2,ijunk
752 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
753 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
756 ! For printing parameters after they are read set the following in the UNRES
759 ! setenv PRINT_PARM YES
761 ! To print parameters in LaTeX format rather than as ASCII tables:
765 call getenv_loc("PRINT_PARM",lancuch)
766 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 call getenv_loc("LATEX",lancuch)
768 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
770 dwa16=2.0d0**(1.0d0/6.0d0)
772 ! Assign virtual-bond length
775 vblinv2=vblinv*vblinv
777 ! Read the virtual-bond parameters, masses, and moments of inertia
778 ! and Stokes' radii of the peptide group and side chains
780 allocate(dsc(ntyp1)) !(ntyp1)
781 allocate(dsc_inv(ntyp1)) !(ntyp1)
782 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
783 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
784 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
785 allocate(nbondterm(ntyp)) !(ntyp)
786 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 allocate(msc(ntyp+1)) !(ntyp+1)
796 allocate(isc(ntyp+1)) !(ntyp+1)
797 allocate(restok(ntyp+1)) !(ntyp+1)
799 read (ibond,*) vbldp0,akp,mp,ip,pstok
802 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
803 dsc(i) = vbldsc0(1,i)
807 dsc_inv(i)=1.0D0/dsc(i)
816 allocate(msc(ntyp+1,5)) !(ntyp+1)
817 allocate(isc(ntyp+1,5)) !(ntyp+1)
818 allocate(restok(ntyp+1,5)) !(ntyp+1)
820 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
822 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
823 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
824 dsc(i) = vbldsc0(1,i)
828 dsc_inv(i)=1.0D0/dsc(i)
832 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
835 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
836 ! dsc(i) = vbldsc0_nucl(1,i)
840 ! dsc_inv(i)=1.0D0/dsc(i)
843 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
844 ! do i=1,ntyp_molec(2)
845 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
846 ! aksc_nucl(j,i),abond0_nucl(j,i),&
847 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
848 ! dsc(i) = vbldsc0(1,i)
852 ! dsc_inv(i)=1.0D0/dsc(i)
857 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
858 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
860 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
862 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
863 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
865 write (iout,'(13x,3f10.5)') &
866 vbldsc0(j,i),aksc(j,i),abond0(j,i)
873 if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5)))
874 if (oldion.eq.1) then
876 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
877 print *,msc(i,5),restok(i,5)
881 read (iion,*) ncatprotparm
882 allocate(catprm(ncatprotparm,4))
884 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
888 allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
892 read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
895 write(*,'(3(5x,a6)11(7x,a6))') "w1 ","w2 ","epslj ","pis1 ", &
896 "sigma0","epsi0 ","chi1 ","chip1 ","sig ","b1 ","b2 ", &
900 write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
901 restyp(i,5),"-",restyp(j,2)
904 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
905 !----------------------------------------------------
906 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
907 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
908 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
909 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
910 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
911 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
921 allocate(liptranene(ntyp))
922 !C reading lipid parameters
923 write (iout,*) "iliptranpar",iliptranpar
925 read(iliptranpar,*) pepliptran
928 read(iliptranpar,*) liptranene(i)
929 print *,liptranene(i)
935 ! Read the parameters of the probability distribution/energy expression
936 ! of the virtual-bond valence angles theta
939 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
940 (bthet(j,i,1,1),j=1,2)
941 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
942 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
943 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
947 athet(1,i,1,-1)=athet(1,i,1,1)
948 athet(2,i,1,-1)=athet(2,i,1,1)
949 bthet(1,i,1,-1)=-bthet(1,i,1,1)
950 bthet(2,i,1,-1)=-bthet(2,i,1,1)
951 athet(1,i,-1,1)=-athet(1,i,1,1)
952 athet(2,i,-1,1)=-athet(2,i,1,1)
953 bthet(1,i,-1,1)=bthet(1,i,1,1)
954 bthet(2,i,-1,1)=bthet(2,i,1,1)
958 athet(1,i,-1,-1)=athet(1,-i,1,1)
959 athet(2,i,-1,-1)=-athet(2,-i,1,1)
960 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
961 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
962 athet(1,i,-1,1)=athet(1,-i,1,1)
963 athet(2,i,-1,1)=-athet(2,-i,1,1)
964 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
965 bthet(2,i,-1,1)=bthet(2,-i,1,1)
966 athet(1,i,1,-1)=-athet(1,-i,1,1)
967 athet(2,i,1,-1)=athet(2,-i,1,1)
968 bthet(1,i,1,-1)=bthet(1,-i,1,1)
969 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
974 polthet(j,i)=polthet(j,-i)
977 gthet(j,i)=gthet(j,-i)
985 'Parameters of the virtual-bond valence angles:'
986 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
987 ' ATHETA0 ',' A1 ',' A2 ',&
990 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
991 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
993 write (iout,'(/a/9x,5a/79(1h-))') &
994 'Parameters of the expression for sigma(theta_c):',&
995 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
996 ' ALPH3 ',' SIGMA0C '
998 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
999 (polthet(j,i),j=0,3),sigc0(i)
1001 write (iout,'(/a/9x,5a/79(1h-))') &
1002 'Parameters of the second gaussian:',&
1003 ' THETA0 ',' SIGMA0 ',' G1 ',&
1006 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
1007 sig0(i),(gthet(j,i),j=1,3)
1010 write (iout,'(a)') &
1011 'Parameters of the virtual-bond valence angles:'
1012 write (iout,'(/a/9x,5a/79(1h-))') &
1013 'Coefficients of expansion',&
1014 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
1015 ' b1*10^1 ',' b2*10^1 '
1017 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
1018 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
1019 (10*bthet(j,i,1,1),j=1,2)
1021 write (iout,'(/a/9x,5a/79(1h-))') &
1022 'Parameters of the expression for sigma(theta_c):',&
1023 ' alpha0 ',' alph1 ',' alph2 ',&
1024 ' alhp3 ',' sigma0c '
1026 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1027 (polthet(j,i),j=0,3),sigc0(i)
1029 write (iout,'(/a/9x,5a/79(1h-))') &
1030 'Parameters of the second gaussian:',&
1031 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1034 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1035 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1041 ! Read the parameters of Utheta determined from ab initio surfaces
1042 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1044 IF (tor_mode.eq.0) THEN
1045 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1046 ntheterm3,nsingle,ndouble
1047 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1049 !----------------------------------------------------
1050 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1051 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1052 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1053 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1054 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1055 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1056 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1057 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1058 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1059 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1060 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1061 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1062 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1063 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1064 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1065 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1066 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1067 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1068 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1069 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1070 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1071 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1072 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1073 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1075 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1077 ithetyp(i)=-ithetyp(-i)
1080 aa0thet(:,:,:,:)=0.0d0
1081 aathet(:,:,:,:,:)=0.0d0
1082 bbthet(:,:,:,:,:,:)=0.0d0
1083 ccthet(:,:,:,:,:,:)=0.0d0
1084 ddthet(:,:,:,:,:,:)=0.0d0
1085 eethet(:,:,:,:,:,:)=0.0d0
1086 ffthet(:,:,:,:,:,:,:)=0.0d0
1087 ggthet(:,:,:,:,:,:,:)=0.0d0
1089 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1091 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1092 ! VAR:1=non-glicyne non-proline 2=proline
1093 ! VAR:negative values for D-aminoacid
1095 do j=-nthetyp,nthetyp
1096 do k=-nthetyp,nthetyp
1097 read (ithep,'(6a)',end=111,err=111) res1
1098 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1099 ! VAR: aa0thet is variable describing the average value of Foureir
1100 ! VAR: expansion series
1101 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1102 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1103 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1104 read (ithep,*,end=111,err=111) &
1105 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1106 read (ithep,*,end=111,err=111) &
1107 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1108 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1109 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1110 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1112 read (ithep,*,end=111,err=111) &
1113 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1114 ffthet(lll,llll,ll,i,j,k,iblock),&
1115 ggthet(llll,lll,ll,i,j,k,iblock),&
1116 ggthet(lll,llll,ll,i,j,k,iblock),&
1117 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1122 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1123 ! coefficients of theta-and-gamma-dependent terms are zero.
1124 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1125 ! RECOMENTDED AFTER VERSION 3.3)
1129 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1130 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1132 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1133 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1136 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1138 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1141 ! AND COMMENT THE LOOPS BELOW
1145 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1146 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1148 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1149 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1152 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1154 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1159 ! Substitution for D aminoacids from symmetry.
1162 do j=-nthetyp,nthetyp
1163 do k=-nthetyp,nthetyp
1164 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1166 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1170 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1171 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1172 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1173 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1179 ffthet(llll,lll,ll,i,j,k,iblock)= &
1180 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1181 ffthet(lll,llll,ll,i,j,k,iblock)= &
1182 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1183 ggthet(llll,lll,ll,i,j,k,iblock)= &
1184 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1185 ggthet(lll,llll,ll,i,j,k,iblock)= &
1186 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1195 ! Control printout of the coefficients of virtual-bond-angle potentials
1198 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1203 write (iout,'(//4a)') &
1204 'Type ',onelett(i),onelett(j),onelett(k)
1205 write (iout,'(//a,10x,a)') " l","a[l]"
1206 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1207 write (iout,'(i2,1pe15.5)') &
1208 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1210 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1211 "b",l,"c",l,"d",l,"e",l
1213 write (iout,'(i2,4(1pe15.5))') m,&
1214 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1215 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1219 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1220 "f+",l,"f-",l,"g+",l,"g-",l
1223 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1224 ffthet(n,m,l,i,j,k,iblock),&
1225 ffthet(m,n,l,i,j,k,iblock),&
1226 ggthet(n,m,l,i,j,k,iblock),&
1227 ggthet(m,n,l,i,j,k,iblock)
1238 !C here will be the apropriate recalibrating for D-aminoacid
1239 read (ithep,*,end=121,err=121) nthetyp
1240 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1241 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1242 do i=-nthetyp+1,nthetyp-1
1243 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1244 do j=0,nbend_kcc_Tb(i)
1245 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1249 write (iout,'(a)') &
1250 "Parameters of the valence-only potentials"
1251 do i=-nthetyp+1,nthetyp-1
1252 write (iout,'(2a)') "Type ",toronelet(i)
1253 do k=0,nbend_kcc_Tb(i)
1254 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1260 write (2,*) "Start reading THETA_PDB",ithep_pdb
1262 ! write (2,*) 'i=',i
1263 read (ithep_pdb,*,err=111,end=111) &
1264 a0thet(i),(athet(j,i,1,1),j=1,2),&
1265 (bthet(j,i,1,1),j=1,2)
1266 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1267 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1268 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1269 sigc0(i)=sigc0(i)**2
1272 athet(1,i,1,-1)=athet(1,i,1,1)
1273 athet(2,i,1,-1)=athet(2,i,1,1)
1274 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1275 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1276 athet(1,i,-1,1)=-athet(1,i,1,1)
1277 athet(2,i,-1,1)=-athet(2,i,1,1)
1278 bthet(1,i,-1,1)=bthet(1,i,1,1)
1279 bthet(2,i,-1,1)=bthet(2,i,1,1)
1282 a0thet(i)=a0thet(-i)
1283 athet(1,i,-1,-1)=athet(1,-i,1,1)
1284 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1285 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1286 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1287 athet(1,i,-1,1)=athet(1,-i,1,1)
1288 athet(2,i,-1,1)=-athet(2,-i,1,1)
1289 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1290 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1291 athet(1,i,1,-1)=-athet(1,-i,1,1)
1292 athet(2,i,1,-1)=athet(2,-i,1,1)
1293 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1294 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1295 theta0(i)=theta0(-i)
1299 polthet(j,i)=polthet(j,-i)
1302 gthet(j,i)=gthet(j,-i)
1305 write (2,*) "End reading THETA_PDB"
1309 !--------------- Reading theta parameters for nucleic acid-------
1310 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1311 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1312 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1313 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1314 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1315 nthetyp_nucl+1,nthetyp_nucl+1))
1316 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1317 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1318 nthetyp_nucl+1,nthetyp_nucl+1))
1319 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1320 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1321 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1322 nthetyp_nucl+1,nthetyp_nucl+1))
1323 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1324 nthetyp_nucl+1,nthetyp_nucl+1))
1325 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1326 nthetyp_nucl+1,nthetyp_nucl+1))
1327 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1328 nthetyp_nucl+1,nthetyp_nucl+1))
1329 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1330 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1331 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1332 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1333 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1334 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1336 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1337 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1339 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1341 aa0thet_nucl(:,:,:)=0.0d0
1342 aathet_nucl(:,:,:,:)=0.0d0
1343 bbthet_nucl(:,:,:,:,:)=0.0d0
1344 ccthet_nucl(:,:,:,:,:)=0.0d0
1345 ddthet_nucl(:,:,:,:,:)=0.0d0
1346 eethet_nucl(:,:,:,:,:)=0.0d0
1347 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1348 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1353 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1354 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1355 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1356 read (ithep_nucl,*,end=111,err=111) &
1357 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1358 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1359 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1360 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1361 read (ithep_nucl,*,end=111,err=111) &
1362 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1363 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1364 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1369 !-------------------------------------------
1370 allocate(nlob(ntyp1)) !(ntyp1)
1371 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1372 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1373 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1380 gaussc(:,:,:,:)=0.0D0
1384 ! Read the parameters of the probability distribution/energy expression
1385 ! of the side chains.
1388 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1392 dsc_inv(i)=1.0D0/dsc(i)
1403 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1404 ((blower(k,l,1),l=1,k),k=1,3)
1405 censc(1,1,-i)=censc(1,1,i)
1406 censc(2,1,-i)=censc(2,1,i)
1407 censc(3,1,-i)=-censc(3,1,i)
1409 read (irotam,*,end=112,err=112) bsc(j,i)
1410 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1411 ((blower(k,l,j),l=1,k),k=1,3)
1412 censc(1,j,-i)=censc(1,j,i)
1413 censc(2,j,-i)=censc(2,j,i)
1414 censc(3,j,-i)=-censc(3,j,i)
1415 ! BSC is amplitude of Gaussian
1422 akl=akl+blower(k,m,j)*blower(l,m,j)
1426 if (((k.eq.3).and.(l.ne.3)) &
1427 .or.((l.eq.3).and.(k.ne.3))) then
1428 gaussc(k,l,j,-i)=-akl
1429 gaussc(l,k,j,-i)=-akl
1431 gaussc(k,l,j,-i)=akl
1432 gaussc(l,k,j,-i)=akl
1441 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1444 if (nlobi.gt.0) then
1446 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1447 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1448 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1449 'log h',(bsc(j,i),j=1,nlobi)
1450 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1451 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1453 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1454 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1457 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1458 write (iout,'(a,f10.4,4(16x,f10.4))') &
1459 'Center ',(bsc(j,i),j=1,nlobi)
1460 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1469 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1470 ! added by Urszula Kozlowska 07/11/2007
1472 !el Maximum number of SC local term fitting function coefficiants
1473 !el integer,parameter :: maxsccoef=65
1475 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1478 read (irotam,*,end=112,err=112)
1480 read (irotam,*,end=112,err=112)
1483 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1487 !---------reading nucleic acid parameters for rotamers-------------------
1488 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1489 do i=1,ntyp_molec(2)
1490 read (irotam_nucl,*,end=112,err=112)
1492 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1498 write (iout,*) "Base rotamer parameters"
1499 do i=1,ntyp_molec(2)
1500 write (iout,'(a)') restyp(i,2)
1501 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1506 ! Read the parameters of the probability distribution/energy expression
1507 ! of the side chains.
1509 write (2,*) "Start reading ROTAM_PDB"
1511 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1515 dsc_inv(i)=1.0D0/dsc(i)
1526 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1527 ((blower(k,l,1),l=1,k),k=1,3)
1529 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1530 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1531 ((blower(k,l,j),l=1,k),k=1,3)
1538 akl=akl+blower(k,m,j)*blower(l,m,j)
1548 write (2,*) "End reading ROTAM_PDB"
1554 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1555 !C interaction energy of the Gly, Ala, and Pro prototypes.
1557 read (ifourier,*) nloctyp
1558 SPLIT_FOURIERTOR = nloctyp.lt.0
1559 nloctyp = iabs(nloctyp)
1560 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1561 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1562 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1563 !C allocate(ctilde(2,2,nres))
1564 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1565 !C allocate(gtb1(2,nres))
1566 !C allocate(gtb2(2,nres))
1567 !C allocate(cc(2,2,nres))
1568 !C allocate(dd(2,2,nres))
1569 !C allocate(ee(2,2,nres))
1570 !C allocate(gtcc(2,2,nres))
1571 !C allocate(gtdd(2,2,nres))
1572 !C allocate(gtee(2,2,nres))
1575 allocate(itype2loc(-ntyp1:ntyp1))
1576 allocate(iloctyp(-nloctyp:nloctyp))
1577 allocate(bnew1(3,2,-nloctyp:nloctyp))
1578 allocate(bnew2(3,2,-nloctyp:nloctyp))
1579 allocate(ccnew(3,2,-nloctyp:nloctyp))
1580 allocate(ddnew(3,2,-nloctyp:nloctyp))
1581 allocate(e0new(3,-nloctyp:nloctyp))
1582 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1583 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1584 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1585 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1586 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1587 allocate(e0newtor(3,-nloctyp:nloctyp))
1588 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1590 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1591 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1592 itype2loc(ntyp1)=nloctyp
1593 iloctyp(nloctyp)=ntyp1
1595 itype2loc(-i)=-itype2loc(i)
1598 allocate(iloctyp(-nloctyp:nloctyp))
1599 allocate(itype2loc(-ntyp1:ntyp1))
1606 iloctyp(-i)=-iloctyp(i)
1608 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1609 !c write (iout,*) "nloctyp",nloctyp,
1610 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1611 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1612 !c write (iout,*) "nloctyp",nloctyp,
1613 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1616 !c write (iout,*) "NEWCORR",i
1617 read (ifourier,*,end=115,err=115)
1620 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1623 !c write (iout,*) "NEWCORR BNEW1"
1624 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1627 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1630 !c write (iout,*) "NEWCORR BNEW2"
1631 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1633 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1634 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1636 !c write (iout,*) "NEWCORR CCNEW"
1637 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1639 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1640 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1642 !c write (iout,*) "NEWCORR DDNEW"
1643 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1647 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1651 !c write (iout,*) "NEWCORR EENEW1"
1652 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1654 read (ifourier,*,end=115,err=115) e0new(ii,i)
1656 !c write (iout,*) (e0new(ii,i),ii=1,3)
1658 !c write (iout,*) "NEWCORR EENEW"
1661 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1662 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1663 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1664 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1669 bnew1(ii,1,-i)= bnew1(ii,1,i)
1670 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1671 bnew2(ii,1,-i)= bnew2(ii,1,i)
1672 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1675 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1676 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1677 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1678 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1679 ccnew(ii,1,-i)=ccnew(ii,1,i)
1680 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1681 ddnew(ii,1,-i)=ddnew(ii,1,i)
1682 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1684 e0new(1,-i)= e0new(1,i)
1685 e0new(2,-i)=-e0new(2,i)
1686 e0new(3,-i)=-e0new(3,i)
1688 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1689 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1690 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1691 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1695 write (iout,'(a)') "Coefficients of the multibody terms"
1696 do i=-nloctyp+1,nloctyp-1
1697 write (iout,*) "Type: ",onelet(iloctyp(i))
1698 write (iout,*) "Coefficients of the expansion of B1"
1700 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1702 write (iout,*) "Coefficients of the expansion of B2"
1704 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1706 write (iout,*) "Coefficients of the expansion of C"
1707 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1708 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1709 write (iout,*) "Coefficients of the expansion of D"
1710 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1711 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1712 write (iout,*) "Coefficients of the expansion of E"
1713 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1716 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1721 IF (SPLIT_FOURIERTOR) THEN
1723 !c write (iout,*) "NEWCORR TOR",i
1724 read (ifourier,*,end=115,err=115)
1727 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1730 !c write (iout,*) "NEWCORR BNEW1 TOR"
1731 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1734 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1737 !c write (iout,*) "NEWCORR BNEW2 TOR"
1738 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1740 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1741 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1743 !c write (iout,*) "NEWCORR CCNEW TOR"
1744 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1746 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1747 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1749 !c write (iout,*) "NEWCORR DDNEW TOR"
1750 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1754 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1758 !c write (iout,*) "NEWCORR EENEW1 TOR"
1759 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1761 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1763 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1765 !c write (iout,*) "NEWCORR EENEW TOR"
1768 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1769 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1770 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1771 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1776 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1777 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1778 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1779 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1782 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1783 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1784 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1785 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1787 e0newtor(1,-i)= e0newtor(1,i)
1788 e0newtor(2,-i)=-e0newtor(2,i)
1789 e0newtor(3,-i)=-e0newtor(3,i)
1791 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1792 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1793 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1794 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1798 write (iout,'(a)') &
1799 "Single-body coefficients of the torsional potentials"
1800 do i=-nloctyp+1,nloctyp-1
1801 write (iout,*) "Type: ",onelet(iloctyp(i))
1802 write (iout,*) "Coefficients of the expansion of B1tor"
1804 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1805 j,(bnew1tor(k,j,i),k=1,3)
1807 write (iout,*) "Coefficients of the expansion of B2tor"
1809 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1810 j,(bnew2tor(k,j,i),k=1,3)
1812 write (iout,*) "Coefficients of the expansion of Ctor"
1813 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1814 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1815 write (iout,*) "Coefficients of the expansion of Dtor"
1816 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1817 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1818 write (iout,*) "Coefficients of the expansion of Etor"
1819 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1822 write (iout,'(1hE,2i1,2f10.5)') &
1823 j,k,(eenewtor(l,j,k,i),l=1,2)
1829 do i=-nloctyp+1,nloctyp-1
1832 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1837 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1841 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1842 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1843 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1844 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1847 ENDIF !SPLIT_FOURIER_TOR
1849 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1850 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1851 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1852 allocate(b(13,-nloctyp-1:nloctyp+1))
1854 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1856 read (ifourier,*,end=115,err=115)
1857 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1859 write (iout,*) 'Type ',onelet(iloctyp(i))
1860 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1874 !c B1tilde(1,i) = b(3)
1875 !c! B1tilde(2,i) =-b(5)
1876 !c! B1tilde(1,-i) =-b(3)
1877 !c! B1tilde(2,-i) =b(5)
1878 !c! b1tilde(1,i)=0.0d0
1879 !c b1tilde(2,i)=0.0d0
1884 !cc B1tilde(1,i) = b(3,i)
1885 !cc B1tilde(2,i) =-b(5,i)
1886 !c B1tilde(1,-i) =-b(3,i)
1887 !c B1tilde(2,-i) =b(5,i)
1888 !cc b1tilde(1,i)=0.0d0
1889 !cc b1tilde(2,i)=0.0d0
1890 !cc B2(1,i) = b(2,i)
1891 !cc B2(2,i) = b(4,i)
1893 !c B2(2,-i) =-b(4,i)
1897 CCold(1,1,i)= b(7,i)
1898 CCold(2,2,i)=-b(7,i)
1899 CCold(2,1,i)= b(9,i)
1900 CCold(1,2,i)= b(9,i)
1901 CCold(1,1,-i)= b(7,i)
1902 CCold(2,2,-i)=-b(7,i)
1903 CCold(2,1,-i)=-b(9,i)
1904 CCold(1,2,-i)=-b(9,i)
1909 !c Ctilde(1,1,i)= CCold(1,1,i)
1910 !c Ctilde(1,2,i)= CCold(1,2,i)
1911 !c Ctilde(2,1,i)=-CCold(2,1,i)
1912 !c Ctilde(2,2,i)=-CCold(2,2,i)
1917 !c Ctilde(1,1,i)= CCold(1,1,i)
1918 !c Ctilde(1,2,i)= CCold(1,2,i)
1919 !c Ctilde(2,1,i)=-CCold(2,1,i)
1920 !c Ctilde(2,2,i)=-CCold(2,2,i)
1922 !c Ctilde(1,1,i)=0.0d0
1923 !c Ctilde(1,2,i)=0.0d0
1924 !c Ctilde(2,1,i)=0.0d0
1925 !c Ctilde(2,2,i)=0.0d0
1926 DDold(1,1,i)= b(6,i)
1927 DDold(2,2,i)=-b(6,i)
1928 DDold(2,1,i)= b(8,i)
1929 DDold(1,2,i)= b(8,i)
1930 DDold(1,1,-i)= b(6,i)
1931 DDold(2,2,-i)=-b(6,i)
1932 DDold(2,1,-i)=-b(8,i)
1933 DDold(1,2,-i)=-b(8,i)
1938 !c Dtilde(1,1,i)= DD(1,1,i)
1939 !c Dtilde(1,2,i)= DD(1,2,i)
1940 !c Dtilde(2,1,i)=-DD(2,1,i)
1941 !c Dtilde(2,2,i)=-DD(2,2,i)
1943 !c Dtilde(1,1,i)=0.0d0
1944 !c Dtilde(1,2,i)=0.0d0
1945 !c Dtilde(2,1,i)=0.0d0
1946 !c Dtilde(2,2,i)=0.0d0
1947 EEold(1,1,i)= b(10,i)+b(11,i)
1948 EEold(2,2,i)=-b(10,i)+b(11,i)
1949 EEold(2,1,i)= b(12,i)-b(13,i)
1950 EEold(1,2,i)= b(12,i)+b(13,i)
1951 EEold(1,1,-i)= b(10,i)+b(11,i)
1952 EEold(2,2,-i)=-b(10,i)+b(11,i)
1953 EEold(2,1,-i)=-b(12,i)+b(13,i)
1954 EEold(1,2,-i)=-b(12,i)-b(13,i)
1955 write(iout,*) "TU DOCHODZE"
1961 !c ee(2,1,i)=ee(1,2,i)
1966 "Coefficients of the cumulants (independent of valence angles)"
1967 do i=-nloctyp+1,nloctyp-1
1968 write (iout,*) 'Type ',onelet(iloctyp(i))
1970 write(iout,'(2f10.5)') B(3,i),B(5,i)
1972 write(iout,'(2f10.5)') B(2,i),B(4,i)
1975 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1979 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1983 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1992 ! Read torsional parameters in old format
1994 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1996 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1997 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1998 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2000 !el from energy module--------
2001 allocate(v1(nterm_old,ntortyp,ntortyp))
2002 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
2003 !el---------------------------
2008 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
2014 write (iout,'(/a/)') 'Torsional constants:'
2017 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
2018 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2024 ! Read torsional parameters
2026 IF (TOR_MODE.eq.0) THEN
2027 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2028 read (itorp,*,end=113,err=113) ntortyp
2029 !el from energy module---------
2030 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2031 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2033 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2034 allocate(vlor2(maxlor,ntortyp,ntortyp))
2035 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2036 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2038 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2039 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2040 !el---------------------------
2043 !el---------------------------
2045 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2047 itortyp(i)=-itortyp(-i)
2049 itortyp(ntyp1)=ntortyp
2050 itortyp(-ntyp1)=-ntortyp
2052 write (iout,*) 'ntortyp',ntortyp
2054 do j=-ntortyp+1,ntortyp-1
2055 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2057 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2058 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2061 do k=1,nterm(i,j,iblock)
2062 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2064 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2065 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2066 v0ij=v0ij+si*v1(k,i,j,iblock)
2068 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2069 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2070 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2072 do k=1,nlor(i,j,iblock)
2073 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2074 vlor2(k,i,j),vlor3(k,i,j)
2075 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2078 v0(-i,-j,iblock)=v0ij
2084 write (iout,'(/a/)') 'Torsional constants:'
2086 do i=-ntortyp,ntortyp
2087 do j=-ntortyp,ntortyp
2088 write (iout,*) 'ityp',i,' jtyp',j
2089 write (iout,*) 'Fourier constants'
2090 do k=1,nterm(i,j,iblock)
2091 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2094 write (iout,*) 'Lorenz constants'
2095 do k=1,nlor(i,j,iblock)
2096 write (iout,'(3(1pe15.5))') &
2097 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2103 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2105 ! 6/23/01 Read parameters for double torsionals
2107 !el from energy module------------
2108 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2109 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2110 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2111 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2112 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2113 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2114 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2115 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2116 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2117 !---------------------------------
2121 do j=-ntortyp+1,ntortyp-1
2122 do k=-ntortyp+1,ntortyp-1
2123 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2124 ! write (iout,*) "OK onelett",
2127 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2128 .or. t3.ne.toronelet(k)) then
2129 write (iout,*) "Error in double torsional parameter file",&
2132 call MPI_Finalize(Ierror)
2134 stop "Error in double torsional parameter file"
2136 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2137 ntermd_2(i,j,k,iblock)
2138 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2139 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2140 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2141 ntermd_1(i,j,k,iblock))
2142 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2143 ntermd_1(i,j,k,iblock))
2144 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2145 ntermd_1(i,j,k,iblock))
2146 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2147 ntermd_1(i,j,k,iblock))
2148 ! Martix of D parameters for one dimesional foureir series
2149 do l=1,ntermd_1(i,j,k,iblock)
2150 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2151 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2152 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2153 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2154 ! write(iout,*) "whcodze" ,
2155 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2157 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2158 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2159 v2s(m,l,i,j,k,iblock),&
2160 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2161 ! Martix of D parameters for two dimesional fourier series
2162 do l=1,ntermd_2(i,j,k,iblock)
2164 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2165 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2166 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2167 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2176 write (iout,*) 'Constants for double torsionals'
2179 do j=-ntortyp+1,ntortyp-1
2180 do k=-ntortyp+1,ntortyp-1
2181 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2182 ' nsingle',ntermd_1(i,j,k,iblock),&
2183 ' ndouble',ntermd_2(i,j,k,iblock)
2185 write (iout,*) 'Single angles:'
2186 do l=1,ntermd_1(i,j,k,iblock)
2187 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2188 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2189 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2190 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2193 write (iout,*) 'Pairs of angles:'
2194 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2195 do l=1,ntermd_2(i,j,k,iblock)
2196 write (iout,'(i5,20f10.5)') &
2197 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2200 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2201 do l=1,ntermd_2(i,j,k,iblock)
2202 write (iout,'(i5,20f10.5)') &
2203 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2204 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2214 itype2loc(i)=itortyp(i)
2218 ELSE IF (TOR_MODE.eq.1) THEN
2220 !C read valence-torsional parameters
2221 read (itorp,*,end=121,err=121) ntortyp
2223 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2225 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2226 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2229 itype2loc(i)=itortyp(i)
2233 itortyp(i)=-itortyp(-i)
2235 do i=-ntortyp+1,ntortyp-1
2236 do j=-ntortyp+1,ntortyp-1
2237 !C first we read the cos and sin gamma parameters
2238 read (itorp,'(13x,a)',end=121,err=121) string
2239 write (iout,*) i,j,string
2240 read (itorp,*,end=121,err=121) &
2241 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2242 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2243 do k=1,nterm_kcc(j,i)
2244 do l=1,nterm_kcc_Tb(j,i)
2245 do ll=1,nterm_kcc_Tb(j,i)
2246 read (itorp,*,end=121,err=121) ii,jj,kk, &
2247 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2255 !c AL 4/8/16: Calculate coefficients from one-body parameters
2257 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2258 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2259 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2260 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2261 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2264 print *,i,itortyp(i)
2265 itortyp(i)=itype2loc(i)
2268 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2270 do i=-ntortyp+1,ntortyp-1
2271 do j=-ntortyp+1,ntortyp-1
2274 do k=1,nterm_kcc_Tb(j,i)
2275 do l=1,nterm_kcc_Tb(j,i)
2276 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2277 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2278 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2279 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2282 do k=1,nterm_kcc_Tb(j,i)
2283 do l=1,nterm_kcc_Tb(j,i)
2285 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2286 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2287 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2288 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2290 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2291 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2292 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2293 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2297 !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
2301 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2306 if (tor_mode.gt.0 .and. lprint) then
2307 !c Print valence-torsional parameters
2308 write (iout,'(a)') &
2309 "Parameters of the valence-torsional potentials"
2310 do i=-ntortyp+1,ntortyp-1
2311 do j=-ntortyp+1,ntortyp-1
2312 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2313 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2314 do k=1,nterm_kcc(j,i)
2315 do l=1,nterm_kcc_Tb(j,i)
2316 do ll=1,nterm_kcc_Tb(j,i)
2317 write (iout,'(3i5,2f15.4)')&
2318 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2326 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2327 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2328 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2329 !el from energy module---------
2330 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2331 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2333 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2334 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2335 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2336 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2338 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2339 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2340 !el---------------------------
2343 !el--------------------
2344 read (itorp_nucl,*,end=113,err=113) &
2345 (itortyp_nucl(i),i=1,ntyp_molec(2))
2346 ! print *,itortyp_nucl(:)
2347 !c write (iout,*) 'ntortyp',ntortyp
2350 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2351 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2354 do k=1,nterm_nucl(i,j)
2355 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2356 v0ij=v0ij+si*v1_nucl(k,i,j)
2359 do k=1,nlor_nucl(i,j)
2360 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2361 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2362 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2368 ! Read of Side-chain backbone correlation parameters
2369 ! Modified 11 May 2012 by Adasko
2372 read (isccor,*,end=119,err=119) nsccortyp
2374 !el from module energy-------------
2375 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2376 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2377 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2378 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2379 !-----------------------------------
2381 !el from module energy-------------
2382 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2384 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2386 isccortyp(i)=-isccortyp(-i)
2388 iscprol=isccortyp(20)
2389 ! write (iout,*) 'ntortyp',ntortyp
2391 !c maxinter is maximum interaction sites
2392 !el from module energy---------
2393 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2394 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2395 -nsccortyp:nsccortyp))
2396 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2397 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2398 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2399 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2400 !-----------------------------------
2404 read (isccor,*,end=119,err=119) &
2405 nterm_sccor(i,j),nlor_sccor(i,j)
2411 nterm_sccor(-i,j)=nterm_sccor(i,j)
2412 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2413 nterm_sccor(i,-j)=nterm_sccor(i,j)
2414 do k=1,nterm_sccor(i,j)
2415 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2417 if (j.eq.iscprol) then
2418 if (i.eq.isccortyp(10)) then
2419 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2420 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2422 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2423 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2424 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2425 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2426 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2427 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2428 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2429 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2432 if (i.eq.isccortyp(10)) then
2433 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2434 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2436 if (j.eq.isccortyp(10)) then
2437 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2438 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2440 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2441 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2442 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2443 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2444 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2445 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2449 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2450 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2451 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2452 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2455 do k=1,nlor_sccor(i,j)
2456 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2457 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2458 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2459 (1+vlor3sccor(k,i,j)**2)
2461 v0sccor(l,i,j)=v0ijsccor
2462 v0sccor(l,-i,j)=v0ijsccor1
2463 v0sccor(l,i,-j)=v0ijsccor2
2464 v0sccor(l,-i,-j)=v0ijsccor3
2470 !el from module energy-------------
2471 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2473 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2474 ! write (iout,*) 'ntortyp',ntortyp
2476 !c maxinter is maximum interaction sites
2477 !el from module energy---------
2478 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2479 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2480 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2481 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2482 !-----------------------------------
2486 read (isccor,*,end=119,err=119) &
2487 nterm_sccor(i,j),nlor_sccor(i,j)
2491 do k=1,nterm_sccor(i,j)
2492 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2494 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2497 do k=1,nlor_sccor(i,j)
2498 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2499 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2500 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2501 (1+vlor3sccor(k,i,j)**2)
2503 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2512 write (iout,'(/a/)') 'Torsional constants:'
2515 write (iout,*) 'ityp',i,' jtyp',j
2516 write (iout,*) 'Fourier constants'
2517 do k=1,nterm_sccor(i,j)
2518 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2520 write (iout,*) 'Lorenz constants'
2521 do k=1,nlor_sccor(i,j)
2522 write (iout,'(3(1pe15.5))') &
2523 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2530 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2531 ! interaction energy of the Gly, Ala, and Pro prototypes.
2534 ! Read electrostatic-interaction parameters
2539 write (iout,'(/a)') 'Electrostatic interaction constants:'
2540 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2541 'IT','JT','APP','BPP','AEL6','AEL3'
2543 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2544 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2545 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2546 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2551 app (i,j)=epp(i,j)*rri*rri
2552 bpp (i,j)=-2.0D0*epp(i,j)*rri
2553 ael6(i,j)=elpp6(i,j)*4.2D0**6
2554 ael3(i,j)=elpp3(i,j)*4.2D0**3
2556 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2562 ! Read side-chain interaction parameters.
2564 !el from module energy - COMMON.INTERACT-------
2565 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2566 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2567 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2569 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2570 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2571 allocate(epslip(ntyp,ntyp))
2579 !--------------------------------
2581 read (isidep,*,end=117,err=117) ipot,expon
2582 if (ipot.lt.1 .or. ipot.gt.5) then
2583 write (iout,'(2a)') 'Error while reading SC interaction',&
2584 'potential file - unknown potential type.'
2586 call MPI_Finalize(Ierror)
2592 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2593 ', exponents are ',expon,2*expon
2594 ! goto (10,20,30,30,40) ipot
2596 !----------------------- LJ potential ---------------------------------
2598 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2599 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2600 (sigma0(i),i=1,ntyp)
2602 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2603 write (iout,'(a/)') 'The epsilon array:'
2604 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2605 write (iout,'(/a)') 'One-body parameters:'
2606 write (iout,'(a,4x,a)') 'residue','sigma'
2607 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2610 !----------------------- LJK potential --------------------------------
2612 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2613 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2614 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2616 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2617 write (iout,'(a/)') 'The epsilon array:'
2618 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2619 write (iout,'(/a)') 'One-body parameters:'
2620 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2621 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2625 !---------------------- GB or BP potential -----------------------------
2628 ! print *,"I AM in SCELE",scelemode
2629 if (scelemode.eq.0) then
2631 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2633 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2634 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2635 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2636 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2638 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2641 ! For the GB potential convert sigma'**2 into chi'
2644 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2648 write (iout,'(/a/)') 'Parameters of the BP potential:'
2649 write (iout,'(a/)') 'The epsilon array:'
2650 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2651 write (iout,'(/a)') 'One-body parameters:'
2652 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2654 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2655 chip(i),alp(i),i=1,ntyp)
2658 ! print *,ntyp,"NTYP"
2659 allocate(icharge(ntyp1))
2660 ! print *,ntyp,icharge(i)
2662 read (isidep,*) (icharge(i),i=1,ntyp)
2663 print *,ntyp,icharge(i)
2664 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2665 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2666 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2667 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2668 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2669 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2670 allocate(epsintab(ntyp,ntyp))
2671 allocate(dtail(2,ntyp,ntyp))
2672 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2673 allocate(wqdip(2,ntyp,ntyp))
2674 allocate(wstate(4,ntyp,ntyp))
2675 allocate(dhead(2,2,ntyp,ntyp))
2676 allocate(nstate(ntyp,ntyp))
2677 allocate(debaykap(ntyp,ntyp))
2679 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2680 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2684 ! write (*,*) "Im in ALAB", i, " ", j
2686 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2687 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2688 chis(i,j),chis(j,i), & !2 w tej linii
2689 nstate(i,j),(wstate(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
2690 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& ! 4 w tej linii
2691 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2692 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2693 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & ! 5 w tej linii
2694 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2695 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2696 IF ((LaTeX).and.(i.gt.24)) then
2697 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2698 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2699 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2700 chis(i,j),chis(j,i) !2 w tej linii
2702 print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i), wqdip(1,i,j)
2707 IF ((LaTeX).and.(i.gt.24)) then
2708 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2709 dhead(1,1,i,j),dhead(2,1,i,j),& ! 2 w tej linii
2710 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2711 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2712 rborn(i,j),rborn(j,i), & ! 3 w tej linii
2713 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2714 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2721 sigma(i,j) = sigma(j,i)
2722 sigmap1(i,j)=sigmap1(j,i)
2723 sigmap2(i,j)=sigmap2(j,i)
2724 sigiso1(i,j)=sigiso1(j,i)
2725 sigiso2(i,j)=sigiso2(j,i)
2726 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2727 nstate(i,j) = nstate(j,i)
2728 dtail(1,i,j) = dtail(1,j,i)
2729 dtail(2,i,j) = dtail(2,j,i)
2731 alphasur(k,i,j) = alphasur(k,j,i)
2732 wstate(k,i,j) = wstate(k,j,i)
2733 alphiso(k,i,j) = alphiso(k,j,i)
2736 dhead(2,1,i,j) = dhead(1,1,j,i)
2737 dhead(2,2,i,j) = dhead(1,2,j,i)
2738 dhead(1,1,i,j) = dhead(2,1,j,i)
2739 dhead(1,2,i,j) = dhead(2,2,j,i)
2741 epshead(i,j) = epshead(j,i)
2742 sig0head(i,j) = sig0head(j,i)
2745 wqdip(k,i,j) = wqdip(k,j,i)
2748 wquad(i,j) = wquad(j,i)
2749 epsintab(i,j) = epsintab(j,i)
2750 debaykap(i,j)=debaykap(j,i)
2751 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2758 !--------------------- GBV potential -----------------------------------
2760 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2761 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2762 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2763 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2765 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2766 write (iout,'(a/)') 'The epsilon array:'
2767 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2768 write (iout,'(/a)') 'One-body parameters:'
2769 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2770 's||/s_|_^2',' chip ',' alph '
2771 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2772 sigii(i),chip(i),alp(i),i=1,ntyp)
2775 write(iout,*)"Wrong ipot"
2781 !-----------------------------------------------------------------------
2782 ! Calculate the "working" parameters of SC interactions.
2784 !el from module energy - COMMON.INTERACT-------
2785 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2786 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2787 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2788 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2789 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2790 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2792 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2798 if (scelemode.eq.0) then
2807 sc_aa_tube_par(:)=0.0d0
2808 sc_bb_tube_par(:)=0.0d0
2810 !--------------------------------
2815 epslip(i,j)=epslip(j,i)
2818 if (scelemode.eq.0) then
2821 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2822 sigma(j,i)=sigma(i,j)
2823 rs0(i,j)=dwa16*sigma(i,j)
2828 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2829 'Working parameters of the SC interactions:',&
2830 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2835 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2837 ! print *,"SIGMA ZLA?",sigma(i,j)
2845 sigeps=dsign(1.0D0,epsij)
2847 aa_aq(i,j)=epsij*rrij*rrij
2848 ! print *,"ADASKO",epsij,rrij,expon
2849 bb_aq(i,j)=-sigeps*epsij*rrij
2850 aa_aq(j,i)=aa_aq(i,j)
2851 bb_aq(j,i)=bb_aq(i,j)
2852 epsijlip=epslip(i,j)
2853 sigeps=dsign(1.0D0,epsijlip)
2854 epsijlip=dabs(epsijlip)
2855 aa_lip(i,j)=epsijlip*rrij*rrij
2856 bb_lip(i,j)=-sigeps*epsijlip*rrij
2857 aa_lip(j,i)=aa_lip(i,j)
2858 bb_lip(j,i)=bb_lip(i,j)
2859 !C write(iout,*) aa_lip
2860 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2861 sigt1sq=sigma0(i)**2
2862 sigt2sq=sigma0(j)**2
2865 ratsig1=sigt2sq/sigt1sq
2866 ratsig2=1.0D0/ratsig1
2867 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2868 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2869 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2873 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2874 sigmaii(i,j)=rsum_max
2875 sigmaii(j,i)=rsum_max
2877 ! sigmaii(i,j)=r0(i,j)
2878 ! sigmaii(j,i)=r0(i,j)
2880 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2881 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2882 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2883 augm(i,j)=epsij*r_augm**(2*expon)
2884 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2891 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2892 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2893 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2898 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2899 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2900 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2901 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2902 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2903 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2904 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2905 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2906 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2907 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2908 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2909 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2910 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2911 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2912 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2913 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2922 read (isidep_nucl,*) ipot_nucl
2923 ! print *,"TU?!",ipot_nucl
2924 if (ipot_nucl.eq.1) then
2925 do i=1,ntyp_molec(2)
2926 do j=i,ntyp_molec(2)
2927 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2928 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2932 do i=1,ntyp_molec(2)
2933 do j=i,ntyp_molec(2)
2934 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2935 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2936 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2940 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2941 do i=1,ntyp_molec(2)
2942 do j=i,ntyp_molec(2)
2943 rrij=sigma_nucl(i,j)
2947 epsij=4*eps_nucl(i,j)
2948 sigeps=dsign(1.0D0,epsij)
2950 aa_nucl(i,j)=epsij*rrij*rrij
2951 bb_nucl(i,j)=-sigeps*epsij*rrij
2952 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2953 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2954 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2955 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2956 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2957 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2960 aa_nucl(i,j)=aa_nucl(j,i)
2961 bb_nucl(i,j)=bb_nucl(j,i)
2962 ael3_nucl(i,j)=ael3_nucl(j,i)
2963 ael6_nucl(i,j)=ael6_nucl(j,i)
2964 ael63_nucl(i,j)=ael63_nucl(j,i)
2965 ael32_nucl(i,j)=ael32_nucl(j,i)
2966 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2967 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2968 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2969 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2970 eps_nucl(i,j)=eps_nucl(j,i)
2971 sigma_nucl(i,j)=sigma_nucl(j,i)
2972 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2976 write(iout,*) "tube param"
2977 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2978 ccavtubpep,dcavtubpep,tubetranenepep
2979 sigmapeptube=sigmapeptube**6
2980 sigeps=dsign(1.0D0,epspeptube)
2981 epspeptube=dabs(epspeptube)
2982 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2983 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2984 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2986 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2987 ccavtub(i),dcavtub(i),tubetranene(i)
2988 sigmasctube=sigmasctube**6
2989 sigeps=dsign(1.0D0,epssctube)
2990 epssctube=dabs(epssctube)
2991 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2992 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2993 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2995 !-----------------READING SC BASE POTENTIALS-----------------------------
2996 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2997 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2998 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2999 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
3000 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
3001 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
3002 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
3003 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
3004 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
3005 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
3006 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
3007 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
3008 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
3009 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
3010 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
3011 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
3014 do i=1,ntyp_molec(1)
3015 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
3016 write (*,*) "Im in ", i, " ", j
3017 read(isidep_scbase,*) &
3018 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3019 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3020 write(*,*) "eps",eps_scbase(i,j)
3021 read(isidep_scbase,*) &
3022 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3023 chis_scbase(i,j,1),chis_scbase(i,j,2)
3024 read(isidep_scbase,*) &
3025 dhead_scbasei(i,j), &
3026 dhead_scbasej(i,j), &
3027 rborn_scbasei(i,j),rborn_scbasej(i,j)
3028 read(isidep_scbase,*) &
3029 (wdipdip_scbase(k,i,j),k=1,3), &
3030 (wqdip_scbase(k,i,j),k=1,2)
3031 read(isidep_scbase,*) &
3032 alphapol_scbase(i,j), &
3033 epsintab_scbase(i,j)
3036 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
3037 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
3039 do i=1,ntyp_molec(1)
3040 do j=1,ntyp_molec(2)-1
3041 epsij=eps_scbase(i,j)
3042 rrij=sigma_scbase(i,j)
3047 sigeps=dsign(1.0D0,epsij)
3049 aa_scbase(i,j)=epsij*rrij*rrij
3050 bb_scbase(i,j)=-sigeps*epsij*rrij
3053 !-----------------READING PEP BASE POTENTIALS-------------------
3054 allocate(eps_pepbase(ntyp_molec(2)))
3055 allocate(sigma_pepbase(ntyp_molec(2)))
3056 allocate(chi_pepbase(ntyp_molec(2),2))
3057 allocate(chipp_pepbase(ntyp_molec(2),2))
3058 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3059 allocate(sigmap1_pepbase(ntyp_molec(2)))
3060 allocate(sigmap2_pepbase(ntyp_molec(2)))
3061 allocate(chis_pepbase(ntyp_molec(2),2))
3062 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3065 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
3066 write (*,*) "Im in ", i, " ", j
3067 read(isidep_pepbase,*) &
3068 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3069 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3070 write(*,*) "eps",eps_pepbase(j)
3071 read(isidep_pepbase,*) &
3072 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3073 chis_pepbase(j,1),chis_pepbase(j,2)
3074 read(isidep_pepbase,*) &
3075 (wdipdip_pepbase(k,j),k=1,3)
3077 allocate(aa_pepbase(ntyp_molec(2)))
3078 allocate(bb_pepbase(ntyp_molec(2)))
3080 do j=1,ntyp_molec(2)-1
3081 epsij=eps_pepbase(j)
3082 rrij=sigma_pepbase(j)
3087 sigeps=dsign(1.0D0,epsij)
3089 aa_pepbase(j)=epsij*rrij*rrij
3090 bb_pepbase(j)=-sigeps*epsij*rrij
3092 !--------------READING SC PHOSPHATE-------------------------------------
3093 allocate(eps_scpho(ntyp_molec(1)))
3094 allocate(sigma_scpho(ntyp_molec(1)))
3095 allocate(chi_scpho(ntyp_molec(1),2))
3096 allocate(chipp_scpho(ntyp_molec(1),2))
3097 allocate(alphasur_scpho(4,ntyp_molec(1)))
3098 allocate(sigmap1_scpho(ntyp_molec(1)))
3099 allocate(sigmap2_scpho(ntyp_molec(1)))
3100 allocate(chis_scpho(ntyp_molec(1),2))
3101 allocate(wqq_scpho(ntyp_molec(1)))
3102 allocate(wqdip_scpho(2,ntyp_molec(1)))
3103 allocate(alphapol_scpho(ntyp_molec(1)))
3104 allocate(epsintab_scpho(ntyp_molec(1)))
3105 allocate(dhead_scphoi(ntyp_molec(1)))
3106 allocate(rborn_scphoi(ntyp_molec(1)))
3107 allocate(rborn_scphoj(ntyp_molec(1)))
3108 allocate(alphi_scpho(ntyp_molec(1)))
3112 do j=1,ntyp_molec(1) ! without U then we will take T for U
3113 write (*,*) "Im in scpho ", i, " ", j
3114 read(isidep_scpho,*) &
3115 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3116 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3117 write(*,*) "eps",eps_scpho(j)
3118 read(isidep_scpho,*) &
3119 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3120 chis_scpho(j,1),chis_scpho(j,2)
3121 read(isidep_scpho,*) &
3122 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3123 read(isidep_scpho,*) &
3124 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3128 allocate(aa_scpho(ntyp_molec(1)))
3129 allocate(bb_scpho(ntyp_molec(1)))
3131 do j=1,ntyp_molec(1)
3138 sigeps=dsign(1.0D0,epsij)
3140 aa_scpho(j)=epsij*rrij*rrij
3141 bb_scpho(j)=-sigeps*epsij*rrij
3145 read(isidep_peppho,*) &
3146 eps_peppho,sigma_peppho
3147 read(isidep_peppho,*) &
3148 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3149 read(isidep_peppho,*) &
3150 (wqdip_peppho(k),k=1,2)
3158 sigeps=dsign(1.0D0,epsij)
3160 aa_peppho=epsij*rrij*rrij
3161 bb_peppho=-sigeps*epsij*rrij
3164 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3169 ! Define the SC-p interaction constants (hard-coded; old style)
3172 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3174 ! aad(i,1)=0.3D0*4.0D0**12
3175 ! Following line for constants currently implemented
3176 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3177 aad(i,1)=1.5D0*4.0D0**12
3178 ! aad(i,1)=0.17D0*5.6D0**12
3180 ! "Soft" SC-p repulsion
3182 ! Following line for constants currently implemented
3183 ! aad(i,1)=0.3D0*4.0D0**6
3184 ! "Hard" SC-p repulsion
3185 bad(i,1)=3.0D0*4.0D0**6
3186 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3195 ! 8/9/01 Read the SC-p interaction constants from file
3198 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3201 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3202 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3203 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3204 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3208 write (iout,*) "Parameters of SC-p interactions:"
3210 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3211 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3216 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3218 do i=1,ntyp_molec(2)
3219 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3221 do i=1,ntyp_molec(2)
3222 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3223 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3225 r0pp=1.12246204830937298142*5.16158
3231 ! Define the constants of the disulfide bridge
3235 ! Old arbitrary potential - commented out.
3240 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3241 ! energy surface of diethyl disulfide.
3242 ! A. Liwo and U. Kozlowska, 11/24/03
3260 allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
3261 allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
3262 allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
3263 allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
3264 allocate(epsintabcat(ntyp,ntyp))
3265 allocate(dtailcat(2,ntyp,ntyp))
3266 allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
3267 allocate(wqdipcat(2,ntyp,ntyp))
3268 allocate(wstatecat(4,ntyp,ntyp))
3269 allocate(dheadcat(2,2,ntyp,ntyp))
3270 allocate(nstatecat(ntyp,ntyp))
3271 allocate(debaykapcat(ntyp,ntyp))
3273 if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
3274 if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
3275 ! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
3276 if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3277 if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3280 if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5)))
3281 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
3282 if (oldion.eq.0) then
3283 if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
3284 allocate(icharge(1:ntyp1))
3285 read(iion,*) (icharge(i),i=1,ntyp)
3290 do i=1,ntyp_molec(5)
3291 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
3292 print *,msc(i,5),restok(i,5)
3296 do j=1,ntyp_molec(5)
3298 ! do j=1,ntyp_molec(5)
3299 ! write (*,*) "Im in ALAB", i, " ", j
3301 epscat(i,j),sigmacat(i,j), &
3302 ! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
3303 chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
3305 (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
3306 ! chiscat(i,j),chiscat(j,i), &
3307 chis1cat(i,j),chis2cat(i,j), &
3309 nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
3310 dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
3311 dtailcat(1,i,j),dtailcat(2,i,j), &
3312 epsheadcat(i,j),sig0headcat(i,j), &
3314 ! rborncat(i,j),rborncat(j,i),&
3315 rborn1cat(i,j),rborn2cat(i,j),&
3316 (wqdipcat(k,i,j),k=1,2), &
3317 alphapolcat(i,j),alphapolcat(j,i), &
3318 (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
3319 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
3321 ! write (iout,*) 'i= ', i, ' j= ', j
3322 ! write (iout,*) 'epsi0= ', epscat(1,j)
3323 ! write (iout,*) 'sigma0= ', sigmacat(1,j)
3324 ! write (iout,*) 'chi(i,j)= ', chicat(1,j)
3325 ! write (iout,*) 'chi(j,i)= ', chicat(j,1)
3326 ! write (iout,*) 'chip(i,j)= ', chippcat(1,j)
3327 ! write (iout,*) 'chip(j,i)= ', chippcat(j,1)
3328 ! write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3329 ! write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3330 ! write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3331 ! write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3332 ! write (iout,*) 'sig1= ', sigmap1cat(1,j)
3333 ! write (iout,*) 'chis(i,j)= ', chiscat(1,j)
3334 ! write (iout,*) 'chis(j,i)= ', chiscat(j,1)
3335 ! write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3336 ! write (iout,*) 'a1= ', rborncat(j,1)
3337 ! write (iout,*) 'a2= ', rborncat(1,j)
3338 ! write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3339 ! write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3340 ! write (iout,*) 'w1= ', wqdipcat(1,1,j)
3341 ! write (iout,*) 'w2= ', wqdipcat(2,1,j)
3345 ! If ((i.eq.1).and.(j.eq.27)) then
3346 ! write (iout,*) 'SEP'
3347 ! Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3348 ! Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3353 allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
3355 do j=1,ntyp_molec(5)
3359 sigeps=dsign(1.0D0,epsij)
3361 aa_aq_cat(i,j)=epsij*rrij*rrij
3362 bb_aq_cat(i,j)=-sigeps*epsij*rrij
3367 do j=1,ntyp_molec(5)
3369 write (iout,*) 'i= ', i, ' j= ', j
3370 write (iout,*) 'epsi0= ', epscat(i,j)
3371 write (iout,*) 'sigma0= ', sigmacat(i,j)
3372 write (iout,*) 'chi1= ', chi1cat(i,j)
3373 write (iout,*) 'chi1= ', chi2cat(i,j)
3374 write (iout,*) 'chip1= ', chipp1cat(1,j)
3375 write (iout,*) 'chip2= ', chipp2cat(1,j)
3376 write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3377 write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3378 write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3379 write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3380 write (iout,*) 'sig1= ', sigmap1cat(1,j)
3381 write (iout,*) 'sig2= ', sigmap2cat(1,j)
3382 write (iout,*) 'chis1= ', chis1cat(1,j)
3383 write (iout,*) 'chis1= ', chis2cat(1,j)
3384 write (iout,*) 'nstatecat(i,j)= ', nstatecat(1,j)
3385 write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,1,j)
3386 write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3387 write (iout,*) 'dhead2= ', dheadcat(1,2,1,j)
3388 write (iout,*) 'a1= ', rborn1cat(i,j)
3389 write (iout,*) 'a2= ', rborn2cat(i,j)
3390 write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3391 write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3392 write (iout,*) 'alphapol2= ', alphapolcat(j,1)
3393 write (iout,*) 'w1= ', wqdipcat(1,i,j)
3394 write (iout,*) 'w2= ', wqdipcat(2,i,j)
3395 write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(1,j)
3398 If ((i.eq.1).and.(j.eq.27)) then
3399 write (iout,*) 'SEP'
3400 Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3401 Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3411 write (iout,'(/a)') "Disulfide bridge parameters:"
3412 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3413 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3414 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3415 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3418 if (shield_mode.gt.0) then
3419 pi=4.0D0*datan(1.0D0)
3420 !C VSolvSphere the volume of solving sphere
3422 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3423 !C there will be no distinction between proline peptide group and normal peptide
3424 !C group in case of shielding parameters
3425 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3426 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3427 write (iout,*) VSolvSphere,VSolvSphere_div
3428 !C long axis of side chain
3430 long_r_sidechain(i)=vbldsc0(1,i)
3431 ! if (scelemode.eq.0) then
3432 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3433 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3435 ! short_r_sidechain(i)=sigma(i,i)
3437 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3444 111 write (iout,*) "Error reading bending energy parameters."
3446 112 write (iout,*) "Error reading rotamer energy parameters."
3448 113 write (iout,*) "Error reading torsional energy parameters."
3450 114 write (iout,*) "Error reading double torsional energy parameters."
3452 115 write (iout,*) &
3453 "Error reading cumulant (multibody energy) parameters."
3455 116 write (iout,*) "Error reading electrostatic energy parameters."
3457 117 write (iout,*) "Error reading side chain interaction parameters."
3459 118 write (iout,*) "Error reading SCp interaction parameters."
3461 119 write (iout,*) "Error reading SCCOR parameters"
3463 121 write (iout,*) "Error in Czybyshev parameters"
3466 call MPI_Finalize(Ierror)
3470 end subroutine parmread
3472 !-----------------------------------------------------------------------------
3474 !-----------------------------------------------------------------------------
3475 subroutine printmat(ldim,m,n,iout,key,a)
3478 character(len=3),dimension(n) :: key
3479 real(kind=8),dimension(ldim,n) :: a
3481 integer :: i,j,k,m,iout,nlim
3485 write (iout,1000) (key(k),k=i,nlim)
3487 1000 format (/5x,8(6x,a3))
3488 1020 format (/80(1h-)/)
3490 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3493 1010 format (a3,2x,8(f9.4))
3495 end subroutine printmat
3496 !-----------------------------------------------------------------------------
3498 !-----------------------------------------------------------------------------
3500 ! Read the PDB file and convert the peptide geometry into virtual-chain
3503 use energy_data, only: itype,istype
3507 ! use control, only: rescode,sugarcode
3508 ! implicit real*8 (a-h,o-z)
3509 ! include 'DIMENSIONS'
3510 ! include 'COMMON.LOCAL'
3511 ! include 'COMMON.VAR'
3512 ! include 'COMMON.CHAIN'
3513 ! include 'COMMON.INTERACT'
3514 ! include 'COMMON.IOUNITS'
3515 ! include 'COMMON.GEO'
3516 ! include 'COMMON.NAMES'
3517 ! include 'COMMON.CONTROL'
3518 ! include 'COMMON.DISTFIT'
3519 ! include 'COMMON.SETUP'
3520 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3522 logical :: lprn=.true.,fail
3523 real(kind=8),dimension(3) :: e1,e2,e3
3524 real(kind=8) :: dcj,efree_temp
3525 character(len=3) :: seq,res,res2
3526 character(len=5) :: atom
3527 character(len=80) :: card
3528 real(kind=8),dimension(3,20) :: sccor
3529 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3530 integer :: isugar,molecprev,firstion
3531 character*1 :: sugar
3533 real(kind=8),dimension(3) :: ccc
3535 integer,dimension(2,maxres/3) :: hfrag_alloc
3536 integer,dimension(4,maxres/3) :: bfrag_alloc
3537 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3538 real(kind=8),dimension(:,:), allocatable :: c_temporary
3539 integer,dimension(:,:) , allocatable :: itype_temporary
3540 integer,dimension(:),allocatable :: istype_temp
3547 ! write (2,*) "UNRES_PDB",unres_pdb
3567 !-----------------------------
3568 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3569 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3570 if(.not. allocated(istype)) allocate(istype(maxres))
3572 read (ipdbin,'(a80)',end=10) card
3573 write (iout,'(a)') card
3574 if (card(:5).eq.'HELIX') then
3577 read(card(22:25),*) hfrag(1,nhfrag)
3578 read(card(34:37),*) hfrag(2,nhfrag)
3580 if (card(:5).eq.'SHEET') then
3583 read(card(24:26),*) bfrag(1,nbfrag)
3584 read(card(35:37),*) bfrag(2,nbfrag)
3585 !rc----------------------------------------
3586 !rc to be corrected !!!
3587 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3588 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3589 !rc----------------------------------------
3591 if (card(:3).eq.'END') then
3593 else if (card(:3).eq.'TER') then
3598 itype(ires_old,molecule)=ntyp1_molec(molecule)
3599 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3600 nres_molec(molecule)=nres_molec(molecule)+2
3602 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3605 dc(j,ires)=sccor(j,iii)
3608 call sccenter(ires,iii,sccor)
3614 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3615 ! Fish out the ATOM cards.
3616 ! write(iout,*) 'card',card(1:20)
3617 ! print *,"ATU ",card(1:6), CARD(3:6)
3619 if (index(card(1:4),'ATOM').gt.0) then
3620 read (card(12:16),*) atom
3621 ! write (iout,*) "! ",atom," !",ires
3622 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3623 read (card(23:26),*) ires
3624 read (card(18:20),'(a3)') res
3625 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3626 ! & " ires_old",ires_old
3627 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3628 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3629 if (ires-ishift+ishift1.ne.ires_old) then
3630 ! Calculate the CM of the preceding residue.
3631 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3633 ! write (iout,*) "Calculating sidechain center iii",iii
3636 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3639 call sccenter(ires_old,iii,sccor)
3643 ! Start new residue.
3644 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3647 else if (ibeg.eq.1) then
3648 write (iout,*) "BEG ires",ires
3650 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3653 nres_molec(molecule)=nres_molec(molecule)+1
3655 ires=ires-ishift+ishift1
3657 ! write (iout,*) "ishift",ishift," ires",ires,&
3658 ! " ires_old",ires_old
3660 else if (ibeg.eq.2) then
3662 ishift=-ires_old+ires-1 !!!!!
3663 ishift1=ishift1-1 !!!!!
3664 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3665 ires=ires-ishift+ishift1
3666 ! print *,ires,ishift,ishift1
3670 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3671 ires=ires-ishift+ishift1
3674 ! print *,'atom',ires,atom
3675 if (res.eq.'ACE' .or. res.eq.'NHE') then
3678 if (atom.eq.'CA '.or.atom.eq.'N ') then
3680 itype(ires,molecule)=rescode(ires,res,0,molecule)
3682 ! nres_molec(molecule)=nres_molec(molecule)+1
3686 itype(ires,molecule)=rescode(ires,res2,0,molecule)
3687 ! nres_molec(molecule)=nres_molec(molecule)+1
3688 read (card(19:19),'(a1)') sugar
3689 isugar=sugarcode(sugar,ires)
3690 ! if (ibeg.eq.1) then
3694 ! print *,"ires=",ires,istype(ires)
3700 ires=ires-ishift+ishift1
3702 ! write (iout,*) "ires_old",ires_old," ires",ires
3703 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3706 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3707 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3708 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3709 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3710 ! print *,ires,ishift,ishift1
3711 ! write (iout,*) "backbone ",atom
3713 write (iout,'(2i3,2x,a,3f8.3)') &
3714 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3717 nres_molec(molecule)=nres_molec(molecule)+1
3719 sccor(j,iii)=c(j,ires)
3721 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3722 atom.eq."C2'" .or. atom.eq."C3'" &
3723 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3724 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3725 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3726 ! print *,ires,ishift,ishift1
3730 ! sccor(j,iii)=c(j,ires)
3733 c(j,ires)=c(j,ires)+ccc(j)/5.0
3735 print *,counter,molecule
3736 if (counter.eq.5) then
3738 nres_molec(molecule)=nres_molec(molecule)+1
3741 ! sccor(j,iii)=c(j,ires)
3745 ! print *, "ATOM",atom(1:3)
3746 else if (atom.eq."C5'") then
3747 read (card(19:19),'(a1)') sugar
3748 isugar=sugarcode(sugar,ires)
3753 ! print *,ires,istype(ires)
3757 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
3758 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3759 nres_molec(molecule)=nres_molec(molecule)+1
3760 print *,"nres_molec(molecule)",nres_molec(molecule),ires
3764 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3766 else if ((atom.eq."C1'").and.unres_pdb) then
3768 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3769 ! write (*,*) card(23:27),ires,itype(ires,1)
3770 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3771 atom.ne.'N' .and. atom.ne.'C' .and. &
3772 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3773 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3774 .and. atom.ne.'P '.and. &
3775 atom(1:1).ne.'H' .and. &
3776 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3778 ! write (iout,*) "sidechain ",atom
3779 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3780 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3781 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3783 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3786 ! print *,"IONS",ions,card(1:6)
3787 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3788 if (firstion.eq.0) then
3792 dc(j,ires)=sccor(j,iii)
3795 call sccenter(ires,iii,sccor)
3798 read (card(12:16),*) atom
3799 ! print *,"HETATOM", atom
3800 read (card(18:20),'(a3)') res
3801 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3802 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3803 .or.(atom(1:2).eq.'K ')) &
3806 if (molecule.ne.5) molecprev=molecule
3808 nres_molec(molecule)=nres_molec(molecule)+1
3809 print *,"HERE",nres_molec(molecule)
3811 itype(ires,molecule)=rescode(ires,res,0,molecule)
3812 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3816 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3817 if (ires.eq.0) return
3818 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3821 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3823 nres_molec(molecule)=nres_molec(molecule)-2
3824 print *,'I have',nres, nres_molec(:)
3826 do k=1,4 ! ions are without dummy
3827 if (nres_molec(k).eq.0) cycle
3829 ! write (iout,*) i,itype(i,1)
3830 ! if (itype(i,1).eq.ntyp1) then
3831 ! write (iout,*) "dummy",i,itype(i,1)
3833 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3834 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3838 if (itype(i,k).eq.ntyp1_molec(k)) then
3839 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3840 if (itype(i+2,k).eq.0) then
3841 ! print *,"masz sieczke"
3843 if (itype(i+2,j).ne.0) then
3845 itype(i+1,j)=ntyp1_molec(j)
3846 nres_molec(k)=nres_molec(k)-1
3847 nres_molec(j)=nres_molec(j)+1
3853 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3854 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3855 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3856 ! if (unres_pdb) then
3857 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3858 ! print *,i,'tu dochodze'
3859 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3867 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3871 dcj=(c(j,i-2)-c(j,i-3))/2.0
3872 if (dcj.eq.0) dcj=1.23591524223
3877 else !itype(i+1,1).eq.ntyp1
3878 ! if (unres_pdb) then
3879 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3880 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3887 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
3888 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
3892 dcj=(c(j,i+3)-c(j,i+2))/2.0
3893 if (dcj.eq.0) dcj=1.23591524223
3898 endif !itype(i+1,1).eq.ntyp1
3899 endif !itype.eq.ntyp1
3903 ! Calculate the CM of the last side chain.
3907 dc(j,ires)=sccor(j,iii)
3910 call sccenter(ires,iii,sccor)
3916 ! print *,"molecule",molecule
3917 if ((itype(nres,1).ne.10)) then
3919 if (molecule.eq.5) molecule=molecprev
3920 itype(nres,molecule)=ntyp1_molec(molecule)
3921 nres_molec(molecule)=nres_molec(molecule)+1
3922 ! if (unres_pdb) then
3923 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3924 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3931 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3935 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3936 c(j,nres)=c(j,nres-1)+dcj
3937 c(j,2*nres)=c(j,nres)
3941 ! print *,'I have',nres, nres_molec(:)
3943 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3945 if (nres.ne.nres0) then
3946 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3948 stop "Error nres value in WHAM input"
3951 !---------------------------------
3952 !el reallocate tables
3955 ! hfrag_alloc(j,i)=hfrag(j,i)
3958 ! bfrag_alloc(j,i)=bfrag(j,i)
3964 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3965 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3966 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3967 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3971 ! hfrag(j,i)=hfrag_alloc(j,i)
3976 ! bfrag(j,i)=bfrag_alloc(j,i)
3979 !el end reallocate tables
3980 !---------------------------------
3988 c(j,2*nres)=c(j,nres)
3991 if (itype(1,1).eq.ntyp1) then
3995 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3996 call refsys(2,3,4,e1,e2,e3,fail)
4003 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4004 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
4008 dcj=(c(j,4)-c(j,3))/2.0
4014 ! First lets assign correct dummy to correct type of chain
4016 if (itype(1,1).eq.ntyp1) then
4017 if (itype(2,1).eq.0) then
4019 if (itype(2,j).ne.0) then
4021 itype(1,j)=ntyp1_molec(j)
4022 nres_molec(1)=nres_molec(1)-1
4023 nres_molec(j)=nres_molec(j)+1
4030 print *,'I have',nres, nres_molec(:)
4032 ! Copy the coordinates to reference coordinates
4038 ! Calculate internal coordinates.
4040 write (iout,'(/a)') &
4041 "Cartesian coordinates of the reference structure"
4042 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4043 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4045 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4046 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4047 (c(j,ires+nres),j=1,3)
4050 ! znamy już nres wiec mozna alokowac tablice
4051 ! Calculate internal coordinates.
4052 if(me.eq.king.or..not.out1file)then
4053 write (iout,'(a)') &
4054 "Backbone and SC coordinates as read from the PDB"
4056 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
4057 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
4058 (c(j,nres+ires),j=1,3)
4061 ! NOW LETS ROCK! SORTING
4062 allocate(c_temporary(3,2*nres))
4063 allocate(itype_temporary(nres,5))
4064 if (.not.allocated(molnum)) allocate(molnum(nres+1))
4065 if (.not.allocated(istype)) write(iout,*) &
4066 "SOMETHING WRONG WITH ISTYTPE"
4067 allocate(istype_temp(nres))
4068 itype_temporary(:,:)=0
4072 if (itype(i,k).ne.0) then
4074 c_temporary(j,seqalingbegin)=c(j,i)
4075 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
4078 itype_temporary(seqalingbegin,k)=itype(i,k)
4079 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
4080 istype_temp(seqalingbegin)=istype(i)
4081 molnum(seqalingbegin)=k
4082 seqalingbegin=seqalingbegin+1
4088 c(j,i)=c_temporary(j,i)
4093 itype(i,k)=itype_temporary(i,k)
4094 istype(i)=istype_temp(i)
4097 if ((itype(1,1).eq.ntyp1).and.itype(2,5).ne.0) then
4098 ! I have only ions now dummy atoms in the system
4100 itype(1,5)=itype(2,5)
4103 itype(i,5)=itype(i+1,5)
4107 nres_molec(1)=nres_molec(1)-1
4109 ! if (itype(1,1).eq.ntyp1) then
4112 ! if (unres_pdb) then
4113 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4114 ! call refsys(2,3,4,e1,e2,e3,fail)
4121 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4125 ! dcj=(c(j,4)-c(j,3))/2.0
4127 ! c(j,nres+1)=c(j,1)
4133 write (iout,'(/a)') &
4134 "Cartesian coordinates of the reference structure after sorting"
4135 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4136 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4138 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4139 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4140 (c(j,ires+nres),j=1,3)
4144 print *,seqalingbegin,nres
4145 if(.not.allocated(vbld)) then
4146 allocate(vbld(2*nres))
4151 if(.not.allocated(vbld_inv)) then
4152 allocate(vbld_inv(2*nres))
4158 if(.not.allocated(theta)) then
4159 allocate(theta(nres+2))
4163 if(.not.allocated(phi)) allocate(phi(nres+2))
4164 if(.not.allocated(alph)) allocate(alph(nres+2))
4165 if(.not.allocated(omeg)) allocate(omeg(nres+2))
4166 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
4167 if(.not.allocated(phiref)) allocate(phiref(nres+2))
4168 if(.not.allocated(costtab)) allocate(costtab(nres))
4169 if(.not.allocated(sinttab)) allocate(sinttab(nres))
4170 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
4171 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
4172 if(.not.allocated(xxref)) allocate(xxref(nres))
4173 if(.not.allocated(yyref)) allocate(yyref(nres))
4174 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
4175 if(.not.allocated(dc_norm)) then
4176 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
4177 allocate(dc_norm(3,0:2*nres+2))
4180 write(iout,*) "before int_from_cart"
4181 call int_from_cart(.true.,.false.)
4182 call sc_loc_geom(.false.)
4183 write(iout,*) "after int_from_cart"
4187 thetaref(i)=theta(i)
4190 write(iout,*) "after thetaref"
4198 dc(j,i)=c(j,i+1)-c(j,i)
4199 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
4204 dc(j,i+nres)=c(j,i+nres)-c(j,i)
4205 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
4207 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
4211 ! Copy the coordinates to reference coordinates
4212 ! Splits to single chain if occurs
4216 ! cref(j,i,cou)=c(j,i)
4222 ! chomo(j,i,k)=c(j,i)
4225 ! write(iout,*) "after chomo"
4227 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
4228 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4229 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4230 !-----------------------------
4234 write (iout,*) "symetr", symetr
4237 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4239 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4242 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4247 cref(j,i,cou)=c(j,i)
4248 cref(j,i+nres,cou)=c(j,i+nres)
4250 chain_rep(j,lll,kkk)=c(j,i)
4251 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4255 write (iout,*) chain_length
4256 if (chain_length.eq.0) chain_length=nres
4258 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
4259 chain_rep(j,chain_length+nres,symetr) &
4260 =chain_rep(j,chain_length+nres,1)
4263 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4265 ! do kkk=1,chain_length
4266 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4270 ! makes copy of chains
4271 write (iout,*) "symetr", symetr
4276 if (symetr.gt.1) then
4283 write(iout,*) (tabperm(i,kkk),kkk=1,4)
4289 write (iout,*) i,icha
4290 do lll=1,chain_length
4292 if (cou.le.nres) then
4294 kupa=mod(lll,chain_length)
4295 iprzes=(kkk-1)*chain_length+lll
4296 if (kupa.eq.0) kupa=chain_length
4297 write (iout,*) "kupa", kupa
4298 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4299 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4306 !-koniec robienia kopii
4309 write (iout,*) "nowa struktura", nperm
4311 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4313 cref(3,i,kkk),cref(1,nres+i,kkk),&
4314 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4316 100 format (//' alpha-carbon coordinates ',&
4317 ' centroid coordinates'/ &
4318 ' ', 6X,'X',11X,'Y',11X,'Z', &
4319 10X,'X',11X,'Y',11X,'Z')
4320 110 format (a,'(',i5,')',6f12.5)
4326 bfrag(i,j)=bfrag(i,j)-ishift
4332 hfrag(i,j)=hfrag(i,j)-ishift
4338 end subroutine readpdb
4339 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4340 !-----------------------------------------------------------------------------
4342 !-----------------------------------------------------------------------------
4343 subroutine read_control
4357 use random, only: random_init
4358 ! implicit real*8 (a-h,o-z)
4359 ! include 'DIMENSIONS'
4361 use prng, only:prng_restart
4363 logical :: OKRandom!, prng_restart
4366 ! include 'COMMON.IOUNITS'
4367 ! include 'COMMON.TIME1'
4368 ! include 'COMMON.THREAD'
4369 ! include 'COMMON.SBRIDGE'
4370 ! include 'COMMON.CONTROL'
4371 ! include 'COMMON.MCM'
4372 ! include 'COMMON.MAP'
4373 ! include 'COMMON.HEADER'
4374 ! include 'COMMON.CSA'
4375 ! include 'COMMON.CHAIN'
4376 ! include 'COMMON.MUCA'
4377 ! include 'COMMON.MD'
4378 ! include 'COMMON.FFIELD'
4379 ! include 'COMMON.INTERACT'
4380 ! include 'COMMON.SETUP'
4381 !el integer :: KDIAG,ICORFL,IXDR
4382 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4383 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4384 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4385 ! character(len=80) :: ucase
4386 character(len=640) :: controlcard
4388 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4390 usampl=.false. ! the default value of usample should be 0
4391 ! write(iout,*) "KURWA2",usampl
4395 read (INP,'(a)') titel
4396 call card_concat(controlcard,.true.)
4397 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4398 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4399 call reada(controlcard,'SEED',seed,0.0D0)
4400 call random_init(seed)
4401 ! Set up the time limit (caution! The time must be input in minutes!)
4402 read_cart=index(controlcard,'READ_CART').gt.0
4403 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4404 call readi(controlcard,'SYM',symetr,1)
4405 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4406 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4407 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4408 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4409 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4410 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4411 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4412 call reada(controlcard,'DRMS',drms,0.1D0)
4413 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
4414 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
4415 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
4416 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
4417 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4418 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4419 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4420 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4421 write (iout,'(a,f10.1)')'DRMS = ',drms
4422 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4423 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4425 call readi(controlcard,'NZ_START',nz_start,0)
4426 call readi(controlcard,'NZ_END',nz_end,0)
4427 call readi(controlcard,'IZ_SC',iz_sc,0)
4428 timlim=60.0D0*timlim
4429 safety = 60.0d0*safety
4432 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4433 !C SHIELD keyword sets if the shielding effect of side-chains is used
4434 !C 0 denots no shielding is used all peptide are equally despite the
4435 !C solvent accesible area
4436 !C 1 the newly introduced function
4437 !C 2 reseved for further possible developement
4438 call readi(controlcard,'SHIELD',shield_mode,0)
4439 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4440 write(iout,*) "shield_mode",shield_mode
4441 call readi(controlcard,'TORMODE',tor_mode,0)
4442 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4443 write(iout,*) "torsional and valence angle mode",tor_mode
4445 !C Varibles set size of box
4446 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4447 protein=index(controlcard,"PROTEIN").gt.0
4448 ions=index(controlcard,"IONS").gt.0
4449 call readi(controlcard,'OLDION',oldion,1)
4450 nucleic=index(controlcard,"NUCLEIC").gt.0
4451 write (iout,*) "with_theta_constr ",with_theta_constr
4452 AFMlog=(index(controlcard,'AFM'))
4453 selfguide=(index(controlcard,'SELFGUIDE'))
4454 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4455 call readi(controlcard,'GENCONSTR',genconstr,0)
4456 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4457 call reada(controlcard,'BOXY',boxysize,100.0d0)
4458 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4459 call readi(controlcard,'TUBEMOD',tubemode,0)
4460 print *,"SCELE",scelemode
4461 call readi(controlcard,"SCELEMODE",scelemode,0)
4462 print *,"SCELE",scelemode
4464 ! elemode = 0 is orignal UNRES electrostatics
4465 ! elemode = 1 is "Momo" potentials in progress
4466 ! elemode = 2 is in development EVALD
4469 write (iout,*) TUBEmode,"TUBEMODE"
4470 if (TUBEmode.gt.0) then
4471 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4472 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4473 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4474 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4475 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4476 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4477 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4478 buftubebot=bordtubebot+tubebufthick
4479 buftubetop=bordtubetop-tubebufthick
4482 ! CUTOFFF ON ELECTROSTATICS
4483 call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
4484 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4485 write(iout,*) "R_CUT_ELE=",r_cut_ele
4486 ! Lipidic parameters
4487 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4488 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4489 if (lipthick.gt.0.0d0) then
4490 bordliptop=(boxzsize+lipthick)/2.0
4491 bordlipbot=bordliptop-lipthick
4492 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4493 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4494 buflipbot=bordlipbot+lipbufthick
4495 bufliptop=bordliptop-lipbufthick
4496 if ((lipbufthick*2.0d0).gt.lipthick) &
4497 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4498 endif !lipthick.gt.0
4499 write(iout,*) "bordliptop=",bordliptop
4500 write(iout,*) "bordlipbot=",bordlipbot
4501 write(iout,*) "bufliptop=",bufliptop
4502 write(iout,*) "buflipbot=",buflipbot
4503 write (iout,*) "SHIELD MODE",shield_mode
4505 !C-------------------------
4506 minim=(index(controlcard,'MINIMIZE').gt.0)
4507 dccart=(index(controlcard,'CART').gt.0)
4508 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4509 overlapsc=.not.overlapsc
4510 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4511 searchsc=.not.searchsc
4512 sideadd=(index(controlcard,'SIDEADD').gt.0)
4513 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4514 outpdb=(index(controlcard,'PDBOUT').gt.0)
4515 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4516 pdbref=(index(controlcard,'PDBREF').gt.0)
4517 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4518 indpdb=index(controlcard,'PDBSTART')
4519 extconf=(index(controlcard,'EXTCONF').gt.0)
4520 call readi(controlcard,'IPRINT',iprint,0)
4521 call readi(controlcard,'MAXGEN',maxgen,10000)
4522 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4523 call readi(controlcard,"KDIAG",kdiag,0)
4524 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4525 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4526 write (iout,*) "RESCALE_MODE",rescale_mode
4527 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4528 if (index(controlcard,'REGULAR').gt.0.0D0) then
4529 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4533 if (index(controlcard,'CHECKGRAD').gt.0) then
4535 if (index(controlcard,'CART').gt.0) then
4537 elseif (index(controlcard,'CARINT').gt.0) then
4542 elseif (index(controlcard,'THREAD').gt.0) then
4544 call readi(controlcard,'THREAD',nthread,0)
4545 if (nthread.gt.0) then
4546 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4549 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4550 stop 'Error termination in Read_Control.'
4552 else if (index(controlcard,'MCMA').gt.0) then
4554 else if (index(controlcard,'MCEE').gt.0) then
4556 else if (index(controlcard,'MULTCONF').gt.0) then
4558 else if (index(controlcard,'MAP').gt.0) then
4560 call readi(controlcard,'MAP',nmap,0)
4561 else if (index(controlcard,'CSA').gt.0) then
4563 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4565 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4568 !fcm else if (index(controlcard,'MCMF').gt.0) then
4570 else if (index(controlcard,'SOFTREG').gt.0) then
4572 else if (index(controlcard,'CHECK_BOND').gt.0) then
4574 else if (index(controlcard,'TEST').gt.0) then
4576 else if (index(controlcard,'MD').gt.0) then
4578 else if (index(controlcard,'RE ').gt.0) then
4582 lmuca=index(controlcard,'MUCA').gt.0
4583 call readi(controlcard,'MUCADYN',mucadyn,0)
4584 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4585 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4587 write (iout,*) 'MUCADYN=',mucadyn
4588 write (iout,*) 'MUCASMOOTH=',muca_smooth
4591 iscode=index(controlcard,'ONE_LETTER')
4592 indphi=index(controlcard,'PHI')
4593 indback=index(controlcard,'BACK')
4594 iranconf=index(controlcard,'RAND_CONF')
4595 i2ndstr=index(controlcard,'USE_SEC_PRED')
4596 gradout=index(controlcard,'GRADOUT').gt.0
4597 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4598 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4599 if (me.eq.king .or. .not.out1file ) &
4600 write (iout,*) "DISTCHAINMAX",distchainmax
4602 if(me.eq.king.or..not.out1file) &
4603 write (iout,'(2a)') diagmeth(kdiag),&
4604 ' routine used to diagonalize matrices.'
4605 if (shield_mode.gt.0) then
4606 pi=4.0D0*datan(1.0D0)
4607 !C VSolvSphere the volume of solving sphere
4609 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4610 !C there will be no distinction between proline peptide group and normal peptide
4611 !C group in case of shielding parameters
4612 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4613 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4614 write (iout,*) VSolvSphere,VSolvSphere_div
4615 !C long axis of side chain
4617 ! long_r_sidechain(i)=vbldsc0(1,i)
4618 ! short_r_sidechain(i)=sigma0(i)
4619 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4626 end subroutine read_control
4627 !-----------------------------------------------------------------------------
4628 subroutine read_REMDpar
4630 ! Read REMD settings
4637 use control_data, only:out1file
4638 ! implicit real*8 (a-h,o-z)
4639 ! include 'DIMENSIONS'
4640 ! include 'COMMON.IOUNITS'
4641 ! include 'COMMON.TIME1'
4642 ! include 'COMMON.MD'
4645 !el include 'COMMON.LANGEVIN'
4647 !el include 'COMMON.LANGEVIN.lang0'
4649 ! include 'COMMON.INTERACT'
4650 ! include 'COMMON.NAMES'
4651 ! include 'COMMON.GEO'
4652 ! include 'COMMON.REMD'
4653 ! include 'COMMON.CONTROL'
4654 ! include 'COMMON.SETUP'
4655 ! character(len=80) :: ucase
4656 character(len=320) :: controlcard
4657 character(len=3200) :: controlcard1
4658 integer :: iremd_m_total
4661 ! real(kind=8) :: var,ene
4663 if(me.eq.king.or..not.out1file) &
4664 write (iout,*) "REMD setup"
4666 call card_concat(controlcard,.true.)
4667 call readi(controlcard,"NREP",nrep,3)
4668 call readi(controlcard,"NSTEX",nstex,1000)
4669 call reada(controlcard,"RETMIN",retmin,10.0d0)
4670 call reada(controlcard,"RETMAX",retmax,1000.0d0)
4671 mremdsync=(index(controlcard,'SYNC').gt.0)
4672 call readi(controlcard,"NSYN",i_sync_step,100)
4673 restart1file=(index(controlcard,'REST1FILE').gt.0)
4674 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
4675 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
4676 if(max_cache_traj_use.gt.max_cache_traj) &
4677 max_cache_traj_use=max_cache_traj
4678 if(me.eq.king.or..not.out1file) then
4679 !d if (traj1file) then
4680 !rc caching is in testing - NTWX is not ignored
4681 !d write (iout,*) "NTWX value is ignored"
4682 !d write (iout,*) " trajectory is stored to one file by master"
4683 !d write (iout,*) " before exchange at NSTEX intervals"
4685 write (iout,*) "NREP= ",nrep
4686 write (iout,*) "NSTEX= ",nstex
4687 write (iout,*) "SYNC= ",mremdsync
4688 write (iout,*) "NSYN= ",i_sync_step
4689 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
4692 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
4693 if (index(controlcard,'TLIST').gt.0) then
4695 call card_concat(controlcard1,.true.)
4696 read(controlcard1,*) (remd_t(i),i=1,nrep)
4697 if(me.eq.king.or..not.out1file) &
4698 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
4701 if (index(controlcard,'MLIST').gt.0) then
4703 call card_concat(controlcard1,.true.)
4704 read(controlcard1,*) (remd_m(i),i=1,nrep)
4705 if(me.eq.king.or..not.out1file) then
4706 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
4709 iremd_m_total=iremd_m_total+remd_m(i)
4711 write (iout,*) 'Total number of replicas ',iremd_m_total
4714 if(me.eq.king.or..not.out1file) &
4715 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4717 end subroutine read_REMDpar
4718 !-----------------------------------------------------------------------------
4719 subroutine read_MDpar
4723 use control_data, only: r_cut,rlamb,out1file
4725 use geometry_data, only: pi
4727 ! implicit real*8 (a-h,o-z)
4728 ! include 'DIMENSIONS'
4729 ! include 'COMMON.IOUNITS'
4730 ! include 'COMMON.TIME1'
4731 ! include 'COMMON.MD'
4734 !el include 'COMMON.LANGEVIN'
4736 !el include 'COMMON.LANGEVIN.lang0'
4738 ! include 'COMMON.INTERACT'
4739 ! include 'COMMON.NAMES'
4740 ! include 'COMMON.GEO'
4741 ! include 'COMMON.SETUP'
4742 ! include 'COMMON.CONTROL'
4743 ! include 'COMMON.SPLITELE'
4744 ! character(len=80) :: ucase
4745 character(len=320) :: controlcard
4750 call card_concat(controlcard,.true.)
4751 call readi(controlcard,"NSTEP",n_timestep,1000000)
4752 call readi(controlcard,"NTWE",ntwe,100)
4753 call readi(controlcard,"NTWX",ntwx,1000)
4754 call reada(controlcard,"DT",d_time,1.0d-1)
4755 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4756 call reada(controlcard,"DAMAX",damax,1.0d1)
4757 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4758 call readi(controlcard,"LANG",lang,0)
4759 RESPA = index(controlcard,"RESPA") .gt. 0
4760 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4761 ntime_split0=ntime_split
4762 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
4763 ntime_split0=ntime_split
4764 call reada(controlcard,"R_CUT",r_cut,2.0d0)
4765 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
4766 rest = index(controlcard,"REST").gt.0
4767 tbf = index(controlcard,"TBF").gt.0
4768 usampl = index(controlcard,"USAMPL").gt.0
4769 ! write(iout,*) "KURWA",usampl
4770 mdpdb = index(controlcard,"MDPDB").gt.0
4771 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4772 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
4773 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
4774 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
4775 if (count_reset_moment.eq.0) count_reset_moment=1000000000
4776 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
4777 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
4778 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
4779 if (count_reset_vel.eq.0) count_reset_vel=1000000000
4780 large = index(controlcard,"LARGE").gt.0
4781 print_compon = index(controlcard,"PRINT_COMPON").gt.0
4782 rattle = index(controlcard,"RATTLE").gt.0
4783 preminim=(index(controlcard,'PREMINIM').gt.0)
4784 write (iout,*) "PREMINIM ",preminim
4785 dccart=(index(controlcard,'CART').gt.0)
4786 if (preminim) call read_minim
4787 ! if performing umbrella sampling, fragments constrained are read from the fragment file
4793 if(me.eq.king.or..not.out1file) then
4795 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
4797 write (iout,'(a)') "The units are:"
4798 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
4799 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
4800 " acceleration: angstrom/(48.9 fs)**2"
4801 write (iout,'(a)') "energy: kcal/mol, temperature: K"
4803 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
4804 write (iout,'(a60,f10.5,a)') &
4805 "Initial time step of numerical integration:",d_time,&
4807 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
4809 write (iout,'(2a,i4,a)') &
4810 "A-MTS algorithm used; initial time step for fast-varying",&
4811 " short-range forces split into",ntime_split," steps."
4812 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
4813 r_cut," lambda",rlamb
4815 write (iout,'(2a,f10.5)') &
4816 "Maximum acceleration threshold to reduce the time step",&
4817 "/increase split number:",damax
4818 write (iout,'(2a,f10.5)') &
4819 "Maximum predicted energy drift to reduce the timestep",&
4820 "/increase split number:",edriftmax
4821 write (iout,'(a60,f10.5)') &
4822 "Maximum velocity threshold to reduce velocities:",dvmax
4823 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4824 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4825 if (rattle) write (iout,'(a60)') &
4826 "Rattle algorithm used to constrain the virtual bonds"
4830 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4831 call reada(controlcard,"RWAT",rwat,1.4d0)
4832 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4833 surfarea=index(controlcard,"SURFAREA").gt.0
4834 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4835 if(me.eq.king.or..not.out1file)then
4836 write (iout,'(/a,$)') "Langevin dynamics calculation"
4838 write (iout,'(a/)') &
4839 " with direct integration of Langevin equations"
4840 else if (lang.eq.2) then
4841 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4842 else if (lang.eq.3) then
4843 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4844 else if (lang.eq.4) then
4845 write (iout,'(a/)') " in overdamped mode"
4847 write (iout,'(//a,i5)') &
4848 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4851 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4852 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4853 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4854 write (iout,'(a60,f10.5)') &
4855 "Scaling factor of the friction forces:",scal_fric
4856 if (surfarea) write (iout,'(2a,i10,a)') &
4857 "Friction coefficients will be scaled by solvent-accessible",&
4858 " surface area every",reset_fricmat," steps."
4860 ! Calculate friction coefficients and bounds of stochastic forces
4861 eta=6*pi*cPoise*etawat
4862 if(me.eq.king.or..not.out1file) &
4863 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4866 do j=1,5 !types of molecules
4867 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4868 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4870 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4871 do j=1,5 !types of molecules
4873 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4874 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4878 if(me.eq.king.or..not.out1file)then
4879 write (iout,'(/2a/)') &
4880 "Radii of site types and friction coefficients and std's of",&
4881 " stochastic forces of fully exposed sites"
4882 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4884 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4885 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4889 if(me.eq.king.or..not.out1file)then
4890 write (iout,'(a)') "Berendsen bath calculation"
4891 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4892 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4894 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4895 count_reset_moment," steps"
4897 write (iout,'(a,i10,a)') &
4898 "Velocities will be reset at random every",count_reset_vel,&
4902 if(me.eq.king.or..not.out1file) &
4903 write (iout,'(a31)') "Microcanonical mode calculation"
4905 if(me.eq.king.or..not.out1file)then
4906 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4908 write(iout,*) "MD running with constraints."
4909 write(iout,*) "Equilibration time ", eq_time, " mtus."
4910 write(iout,*) "Constraining ", nfrag," fragments."
4911 write(iout,*) "Length of each fragment, weight and q0:"
4913 write (iout,*) "Set of restraints #",iset
4915 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4916 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4918 write(iout,*) "constraints between ", npair, "fragments."
4919 write(iout,*) "constraint pairs, weights and q0:"
4921 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4922 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4924 write(iout,*) "angle constraints within ", nfrag_back,&
4925 "backbone fragments."
4926 write(iout,*) "fragment, weights:"
4928 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4929 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4930 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4933 iset=mod(kolor,nset)+1
4936 if(me.eq.king.or..not.out1file) &
4937 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4939 end subroutine read_MDpar
4940 !-----------------------------------------------------------------------------
4944 ! implicit real*8 (a-h,o-z)
4945 ! include 'DIMENSIONS'
4946 ! include 'COMMON.MAP'
4947 ! include 'COMMON.IOUNITS'
4948 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4949 character(len=80) :: mapcard !,ucase
4952 ! real(kind=8) :: var,ene
4955 read (inp,'(a)') mapcard
4956 mapcard=ucase(mapcard)
4957 if (index(mapcard,'PHI').gt.0) then
4959 else if (index(mapcard,'THE').gt.0) then
4961 else if (index(mapcard,'ALP').gt.0) then
4963 else if (index(mapcard,'OME').gt.0) then
4966 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4967 stop 'Error - illegal variable spec in MAP card.'
4969 call readi (mapcard,'RES1',res1(imap),0)
4970 call readi (mapcard,'RES2',res2(imap),0)
4971 if (res1(imap).eq.0) then
4972 res1(imap)=res2(imap)
4973 else if (res2(imap).eq.0) then
4974 res2(imap)=res1(imap)
4976 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4977 write (iout,'(a)') &
4978 'Error - illegal definition of variable group in MAP.'
4979 stop 'Error - illegal definition of variable group in MAP.'
4981 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4982 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4983 call readi(mapcard,'NSTEP',nstep(imap),0)
4984 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4985 write (iout,'(a)') &
4986 'Illegal boundary and/or step size specification in MAP.'
4987 stop 'Illegal boundary and/or step size specification in MAP.'
4991 end subroutine map_read
4992 !-----------------------------------------------------------------------------
4995 use control_data, only: vdisulf
4997 ! implicit real*8 (a-h,o-z)
4998 ! include 'DIMENSIONS'
4999 ! include 'COMMON.IOUNITS'
5000 ! include 'COMMON.GEO'
5001 ! include 'COMMON.CSA'
5002 ! include 'COMMON.BANK'
5003 ! include 'COMMON.CONTROL'
5004 ! character(len=80) :: ucase
5005 character(len=620) :: mcmcard
5007 ! integer :: ntf,ik,iw_pdb
5008 ! real(kind=8) :: var,ene
5010 call card_concat(mcmcard,.true.)
5012 call readi(mcmcard,'NCONF',nconf,50)
5013 call readi(mcmcard,'NADD',nadd,0)
5014 call readi(mcmcard,'JSTART',jstart,1)
5015 call readi(mcmcard,'JEND',jend,1)
5016 call readi(mcmcard,'NSTMAX',nstmax,500000)
5017 call readi(mcmcard,'N0',n0,1)
5018 call readi(mcmcard,'N1',n1,6)
5019 call readi(mcmcard,'N2',n2,4)
5020 call readi(mcmcard,'N3',n3,0)
5021 call readi(mcmcard,'N4',n4,0)
5022 call readi(mcmcard,'N5',n5,0)
5023 call readi(mcmcard,'N6',n6,10)
5024 call readi(mcmcard,'N7',n7,0)
5025 call readi(mcmcard,'N8',n8,0)
5026 call readi(mcmcard,'N9',n9,0)
5027 call readi(mcmcard,'N14',n14,0)
5028 call readi(mcmcard,'N15',n15,0)
5029 call readi(mcmcard,'N16',n16,0)
5030 call readi(mcmcard,'N17',n17,0)
5031 call readi(mcmcard,'N18',n18,0)
5033 vdisulf=(index(mcmcard,'DYNSS').gt.0)
5035 call readi(mcmcard,'NDIFF',ndiff,2)
5036 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
5037 call readi(mcmcard,'IS1',is1,1)
5038 call readi(mcmcard,'IS2',is2,8)
5039 call readi(mcmcard,'NRAN0',nran0,4)
5040 call readi(mcmcard,'NRAN1',nran1,2)
5041 call readi(mcmcard,'IRR',irr,1)
5042 call readi(mcmcard,'NSEED',nseed,20)
5043 call readi(mcmcard,'NTOTAL',ntotal,10000)
5044 call reada(mcmcard,'CUT1',cut1,2.0d0)
5045 call reada(mcmcard,'CUT2',cut2,5.0d0)
5046 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
5047 call readi(mcmcard,'ICMAX',icmax,3)
5048 call readi(mcmcard,'IRESTART',irestart,0)
5049 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
5052 call reada(mcmcard,'DELE',dele,20.0d0)
5053 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
5054 call readi(mcmcard,'IREF',iref,0)
5055 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
5056 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
5057 call readi(mcmcard,'NCONF_IN',nconf_in,0)
5058 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
5059 write (iout,*) "NCONF_IN",nconf_in
5061 end subroutine csaread
5062 !-----------------------------------------------------------------------------
5066 use control_data, only: MaxMoveType
5069 ! implicit real*8 (a-h,o-z)
5070 ! include 'DIMENSIONS'
5071 ! include 'COMMON.MCM'
5072 ! include 'COMMON.MCE'
5073 ! include 'COMMON.IOUNITS'
5074 ! character(len=80) :: ucase
5075 character(len=320) :: mcmcard
5078 ! real(kind=8) :: var,ene
5080 call card_concat(mcmcard,.true.)
5081 call readi(mcmcard,'MAXACC',maxacc,100)
5082 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
5083 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
5084 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
5085 call readi(mcmcard,'MAXREPM',maxrepm,200)
5086 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
5087 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
5088 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
5089 call reada(mcmcard,'E_UP',e_up,5.0D0)
5090 call reada(mcmcard,'DELTE',delte,0.1D0)
5091 call readi(mcmcard,'NSWEEP',nsweep,5)
5092 call readi(mcmcard,'NSTEPH',nsteph,0)
5093 call readi(mcmcard,'NSTEPC',nstepc,0)
5094 call reada(mcmcard,'TMIN',tmin,298.0D0)
5095 call reada(mcmcard,'TMAX',tmax,298.0D0)
5096 call readi(mcmcard,'NWINDOW',nwindow,0)
5097 call readi(mcmcard,'PRINT_MC',print_mc,0)
5098 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
5099 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
5100 ent_read=(index(mcmcard,'ENT_READ').gt.0)
5101 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
5102 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
5103 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
5104 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
5105 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
5106 if (nwindow.gt.0) then
5107 allocate(winstart(nwindow)) !!el (maxres)
5108 allocate(winend(nwindow)) !!el
5109 allocate(winlen(nwindow)) !!el
5110 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
5112 winlen(i)=winend(i)-winstart(i)+1
5115 if (tmax.lt.tmin) tmax=tmin
5116 if (tmax.eq.tmin) then
5120 if (nstepc.gt.0 .and. nsteph.gt.0) then
5121 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
5122 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
5124 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
5125 ! Probabilities of different move types
5126 sumpro_type(0)=0.0D0
5127 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
5128 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
5129 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
5130 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
5131 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
5132 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
5133 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
5135 print *,'i',i,' sumprotype',sumpro_type(i)
5136 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
5137 print *,'i',i,' sumprotype',sumpro_type(i)
5140 end subroutine mcmread
5141 !-----------------------------------------------------------------------------
5142 subroutine read_minim
5145 ! implicit real*8 (a-h,o-z)
5146 ! include 'DIMENSIONS'
5147 ! include 'COMMON.MINIM'
5148 ! include 'COMMON.IOUNITS'
5149 ! character(len=80) :: ucase
5150 character(len=320) :: minimcard
5152 ! integer :: ntf,ik,iw_pdb
5153 ! real(kind=8) :: var,ene
5155 call card_concat(minimcard,.true.)
5156 call readi(minimcard,'MAXMIN',maxmin,2000)
5157 call readi(minimcard,'MAXFUN',maxfun,5000)
5158 call readi(minimcard,'MINMIN',minmin,maxmin)
5159 call readi(minimcard,'MINFUN',minfun,maxmin)
5160 call reada(minimcard,'TOLF',tolf,1.0D-2)
5161 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
5162 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
5163 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
5164 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
5165 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
5166 'Options in energy minimization:'
5167 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
5168 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
5169 'MinMin:',MinMin,' MinFun:',MinFun,&
5170 ' TolF:',TolF,' RTolF:',RTolF
5172 end subroutine read_minim
5173 !-----------------------------------------------------------------------------
5174 subroutine openunits
5176 use MD_data, only: usampl
5179 use control_data, only:out1file
5180 use control, only: getenv_loc
5181 ! implicit real*8 (a-h,o-z)
5182 ! include 'DIMENSIONS'
5185 character(len=16) :: form,nodename
5186 integer :: nodelen,ierror,npos
5188 ! include 'COMMON.SETUP'
5189 ! include 'COMMON.IOUNITS'
5190 ! include 'COMMON.MD'
5191 ! include 'COMMON.CONTROL'
5192 integer :: lenpre,lenpot,lentmp !,ilen
5194 character(len=3) :: out1file_text !,ucase
5195 character(len=3) :: ll
5198 ! integer :: ntf,ik,iw_pdb
5199 ! real(kind=8) :: var,ene
5201 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
5202 call getenv_loc("PREFIX",prefix)
5204 call getenv_loc("POT",pot)
5205 call getenv_loc("DIRTMP",tmpdir)
5206 call getenv_loc("CURDIR",curdir)
5207 call getenv_loc("OUT1FILE",out1file_text)
5208 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
5209 out1file_text=ucase(out1file_text)
5210 if (out1file_text(1:1).eq."Y") then
5213 out1file=fg_rank.gt.0
5218 if (lentmp.gt.0) then
5219 write (*,'(80(1h!))')
5220 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
5221 write (*,'(80(1h!))')
5222 write (*,*)"All output files will be on node /tmp directory."
5224 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
5225 if (me.eq.king) then
5226 write (*,*) "The master node is ",nodename
5227 else if (fg_rank.eq.0) then
5228 write (*,*) "I am the CG slave node ",nodename
5230 write (*,*) "I am the FG slave node ",nodename
5233 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
5234 lenpre = lentmp+lenpre+1
5236 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
5237 ! Get the names and open the input files
5238 #if defined(WINIFL) || defined(WINPGI)
5239 open(1,file=pref_orig(:ilen(pref_orig))// &
5240 '.inp',status='old',readonly,shared)
5241 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5242 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5243 ! Get parameter filenames and open the parameter files.
5244 call getenv_loc('BONDPAR',bondname)
5245 open (ibond,file=bondname,status='old',readonly,shared)
5246 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5247 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5248 call getenv_loc('THETPAR',thetname)
5249 open (ithep,file=thetname,status='old',readonly,shared)
5250 call getenv_loc('ROTPAR',rotname)
5251 open (irotam,file=rotname,status='old',readonly,shared)
5252 call getenv_loc('TORPAR',torname)
5253 open (itorp,file=torname,status='old',readonly,shared)
5254 call getenv_loc('TORDPAR',tordname)
5255 open (itordp,file=tordname,status='old',readonly,shared)
5256 call getenv_loc('FOURIER',fouriername)
5257 open (ifourier,file=fouriername,status='old',readonly,shared)
5258 call getenv_loc('ELEPAR',elename)
5259 open (ielep,file=elename,status='old',readonly,shared)
5260 call getenv_loc('SIDEPAR',sidename)
5261 open (isidep,file=sidename,status='old',readonly,shared)
5263 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5264 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5265 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5266 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5267 call getenv_loc('TORPAR_NUCL',torname_nucl)
5268 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5269 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5270 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5271 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5272 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5275 #elif (defined CRAY) || (defined AIX)
5276 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5278 ! print *,"Processor",myrank," opened file 1"
5279 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5280 ! print *,"Processor",myrank," opened file 9"
5281 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5282 ! Get parameter filenames and open the parameter files.
5283 call getenv_loc('BONDPAR',bondname)
5284 open (ibond,file=bondname,status='old',action='read')
5285 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5286 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5288 ! print *,"Processor",myrank," opened file IBOND"
5289 call getenv_loc('THETPAR',thetname)
5290 open (ithep,file=thetname,status='old',action='read')
5291 ! print *,"Processor",myrank," opened file ITHEP"
5292 call getenv_loc('ROTPAR',rotname)
5293 open (irotam,file=rotname,status='old',action='read')
5294 ! print *,"Processor",myrank," opened file IROTAM"
5295 call getenv_loc('TORPAR',torname)
5296 open (itorp,file=torname,status='old',action='read')
5297 ! print *,"Processor",myrank," opened file ITORP"
5298 call getenv_loc('TORDPAR',tordname)
5299 open (itordp,file=tordname,status='old',action='read')
5300 ! print *,"Processor",myrank," opened file ITORDP"
5301 call getenv_loc('SCCORPAR',sccorname)
5302 open (isccor,file=sccorname,status='old',action='read')
5303 ! print *,"Processor",myrank," opened file ISCCOR"
5304 call getenv_loc('FOURIER',fouriername)
5305 open (ifourier,file=fouriername,status='old',action='read')
5306 ! print *,"Processor",myrank," opened file IFOURIER"
5307 call getenv_loc('ELEPAR',elename)
5308 open (ielep,file=elename,status='old',action='read')
5309 ! print *,"Processor",myrank," opened file IELEP"
5310 call getenv_loc('SIDEPAR',sidename)
5311 open (isidep,file=sidename,status='old',action='read')
5313 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5314 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5315 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5316 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5317 call getenv_loc('TORPAR_NUCL',torname_nucl)
5318 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5319 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5320 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5321 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5322 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5324 call getenv_loc('LIPTRANPAR',liptranname)
5325 open (iliptranpar,file=liptranname,status='old',action='read')
5326 call getenv_loc('TUBEPAR',tubename)
5327 open (itube,file=tubename,status='old',action='read')
5328 call getenv_loc('IONPAR',ionname)
5329 open (iion,file=ionname,status='old',action='read')
5331 ! print *,"Processor",myrank," opened file ISIDEP"
5332 ! print *,"Processor",myrank," opened parameter files"
5334 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5335 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5336 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5337 ! Get parameter filenames and open the parameter files.
5338 call getenv_loc('BONDPAR',bondname)
5339 open (ibond,file=bondname,status='old')
5340 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5341 open (ibond_nucl,file=bondname_nucl,status='old')
5343 call getenv_loc('THETPAR',thetname)
5344 open (ithep,file=thetname,status='old')
5345 call getenv_loc('ROTPAR',rotname)
5346 open (irotam,file=rotname,status='old')
5347 call getenv_loc('TORPAR',torname)
5348 open (itorp,file=torname,status='old')
5349 call getenv_loc('TORDPAR',tordname)
5350 open (itordp,file=tordname,status='old')
5351 call getenv_loc('SCCORPAR',sccorname)
5352 open (isccor,file=sccorname,status='old')
5353 call getenv_loc('FOURIER',fouriername)
5354 open (ifourier,file=fouriername,status='old')
5355 call getenv_loc('ELEPAR',elename)
5356 open (ielep,file=elename,status='old')
5357 call getenv_loc('SIDEPAR',sidename)
5358 open (isidep,file=sidename,status='old')
5360 open (ithep_nucl,file=thetname_nucl,status='old')
5361 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5362 open (irotam_nucl,file=rotname_nucl,status='old')
5363 call getenv_loc('TORPAR_NUCL',torname_nucl)
5364 open (itorp_nucl,file=torname_nucl,status='old')
5365 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5366 open (itordp_nucl,file=tordname_nucl,status='old')
5367 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5368 open (isidep_nucl,file=sidename_nucl,status='old')
5370 call getenv_loc('LIPTRANPAR',liptranname)
5371 open (iliptranpar,file=liptranname,status='old')
5372 call getenv_loc('TUBEPAR',tubename)
5373 open (itube,file=tubename,status='old')
5374 call getenv_loc('IONPAR',ionname)
5375 open (iion,file=ionname,status='old')
5376 call getenv_loc('IONPAR_NUCL',ionnuclname)
5377 open (iionnucl,file=ionnuclname,status='old')
5379 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5381 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5382 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5383 ! Get parameter filenames and open the parameter files.
5384 call getenv_loc('BONDPAR',bondname)
5385 open (ibond,file=bondname,status='old',action='read')
5386 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5387 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5388 call getenv_loc('THETPAR',thetname)
5389 open (ithep,file=thetname,status='old',action='read')
5390 call getenv_loc('ROTPAR',rotname)
5391 open (irotam,file=rotname,status='old',action='read')
5392 call getenv_loc('TORPAR',torname)
5393 open (itorp,file=torname,status='old',action='read')
5394 call getenv_loc('TORDPAR',tordname)
5395 open (itordp,file=tordname,status='old',action='read')
5396 call getenv_loc('SCCORPAR',sccorname)
5397 open (isccor,file=sccorname,status='old',action='read')
5399 call getenv_loc('THETPARPDB',thetname_pdb)
5400 print *,"thetname_pdb ",thetname_pdb
5401 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5402 print *,ithep_pdb," opened"
5404 call getenv_loc('FOURIER',fouriername)
5405 open (ifourier,file=fouriername,status='old',readonly)
5406 call getenv_loc('ELEPAR',elename)
5407 open (ielep,file=elename,status='old',readonly)
5408 call getenv_loc('SIDEPAR',sidename)
5409 open (isidep,file=sidename,status='old',readonly)
5411 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5412 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5413 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5414 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5415 call getenv_loc('TORPAR_NUCL',torname_nucl)
5416 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5417 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5418 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5419 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5420 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5421 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5422 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5423 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5424 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5425 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5426 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5427 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5428 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5431 call getenv_loc('LIPTRANPAR',liptranname)
5432 open (iliptranpar,file=liptranname,status='old',action='read')
5433 call getenv_loc('TUBEPAR',tubename)
5434 open (itube,file=tubename,status='old',action='read')
5435 call getenv_loc('IONPAR',ionname)
5436 open (iion,file=ionname,status='old',action='read')
5437 call getenv_loc('IONPAR_NUCL',ionnuclname)
5438 open (iionnucl,file=ionnuclname,status='old',action='read')
5441 call getenv_loc('ROTPARPDB',rotname_pdb)
5442 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5445 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5446 #if defined(WINIFL) || defined(WINPGI)
5447 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5448 #elif (defined CRAY) || (defined AIX)
5449 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5451 open (iscpp_nucl,file=scpname_nucl,status='old')
5453 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5458 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5459 ! Use -DOLDSCP to use hard-coded constants instead.
5461 call getenv_loc('SCPPAR',scpname)
5462 #if defined(WINIFL) || defined(WINPGI)
5463 open (iscpp,file=scpname,status='old',readonly,shared)
5464 #elif (defined CRAY) || (defined AIX)
5465 open (iscpp,file=scpname,status='old',action='read')
5467 open (iscpp,file=scpname,status='old')
5469 open (iscpp,file=scpname,status='old',action='read')
5472 call getenv_loc('PATTERN',patname)
5473 #if defined(WINIFL) || defined(WINPGI)
5474 open (icbase,file=patname,status='old',readonly,shared)
5475 #elif (defined CRAY) || (defined AIX)
5476 open (icbase,file=patname,status='old',action='read')
5478 open (icbase,file=patname,status='old')
5480 open (icbase,file=patname,status='old',action='read')
5483 ! Open output file only for CG processes
5484 ! print *,"Processor",myrank," fg_rank",fg_rank
5485 if (fg_rank.eq.0) then
5487 if (nodes.eq.1) then
5490 npos = dlog10(dfloat(nodes-1))+1
5492 if (npos.lt.3) npos=3
5493 write (liczba,'(i1)') npos
5494 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5496 write (liczba,form) me
5497 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5498 liczba(:ilen(liczba))
5499 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5501 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5503 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5504 liczba(:ilen(liczba))//'.mol2'
5505 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5506 liczba(:ilen(liczba))//'.stat'
5508 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5509 //liczba(:ilen(liczba))//'.stat')
5510 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5513 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5514 liczba(:ilen(liczba))//'.const'
5519 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5520 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5521 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5522 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5523 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5525 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5527 rest2name=prefix(:ilen(prefix))//'.rst'
5529 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5532 #if defined(AIX) || defined(PGI)
5533 if (me.eq.king .or. .not. out1file) &
5534 open(iout,file=outname,status='unknown')
5536 if (fg_rank.gt.0) then
5537 write (liczba,'(i3.3)') myrank/nfgtasks
5538 write (ll,'(bz,i3.3)') fg_rank
5539 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5544 open(igeom,file=intname,status='unknown',position='append')
5545 open(ipdb,file=pdbname,status='unknown')
5546 open(imol2,file=mol2name,status='unknown')
5547 open(istat,file=statname,status='unknown',position='append')
5549 !1out open(iout,file=outname,status='unknown')
5552 if (me.eq.king .or. .not.out1file) &
5553 open(iout,file=outname,status='unknown')
5555 if (fg_rank.gt.0) then
5556 write (liczba,'(i3.3)') myrank/nfgtasks
5557 write (ll,'(bz,i3.3)') fg_rank
5558 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5563 open(igeom,file=intname,status='unknown',access='append')
5564 open(ipdb,file=pdbname,status='unknown')
5565 open(imol2,file=mol2name,status='unknown')
5566 open(istat,file=statname,status='unknown',access='append')
5568 !1out open(iout,file=outname,status='unknown')
5571 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5572 csa_seed=prefix(:lenpre)//'.CSA.seed'
5573 csa_history=prefix(:lenpre)//'.CSA.history'
5574 csa_bank=prefix(:lenpre)//'.CSA.bank'
5575 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5576 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5577 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5578 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5579 csa_int=prefix(:lenpre)//'.int'
5580 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5581 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5582 csa_in=prefix(:lenpre)//'.CSA.in'
5583 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5586 write (iout,'(80(1h-))')
5587 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5588 write (iout,'(80(1h-))')
5589 write (iout,*) "Input file : ",&
5590 pref_orig(:ilen(pref_orig))//'.inp'
5591 write (iout,*) "Output file : ",&
5592 outname(:ilen(outname))
5594 write (iout,*) "Sidechain potential file : ",&
5595 sidename(:ilen(sidename))
5597 write (iout,*) "SCp potential file : ",&
5598 scpname(:ilen(scpname))
5600 write (iout,*) "Electrostatic potential file : ",&
5601 elename(:ilen(elename))
5602 write (iout,*) "Cumulant coefficient file : ",&
5603 fouriername(:ilen(fouriername))
5604 write (iout,*) "Torsional parameter file : ",&
5605 torname(:ilen(torname))
5606 write (iout,*) "Double torsional parameter file : ",&
5607 tordname(:ilen(tordname))
5608 write (iout,*) "SCCOR parameter file : ",&
5609 sccorname(:ilen(sccorname))
5610 write (iout,*) "Bond & inertia constant file : ",&
5611 bondname(:ilen(bondname))
5612 write (iout,*) "Bending parameter file : ",&
5613 thetname(:ilen(thetname))
5614 write (iout,*) "Rotamer parameter file : ",&
5615 rotname(:ilen(rotname))
5618 write (iout,*) "Thetpdb parameter file : ",&
5619 thetname_pdb(:ilen(thetname_pdb))
5622 write (iout,*) "Threading database : ",&
5623 patname(:ilen(patname))
5625 write (iout,*)" DIRTMP : ",&
5627 write (iout,'(80(1h-))')
5630 end subroutine openunits
5631 !-----------------------------------------------------------------------------
5634 use geometry_data, only: nres,dc
5636 ! implicit real*8 (a-h,o-z)
5637 ! include 'DIMENSIONS'
5638 ! include 'COMMON.CHAIN'
5639 ! include 'COMMON.IOUNITS'
5640 ! include 'COMMON.MD'
5643 ! real(kind=8) :: var,ene
5645 open(irest2,file=rest2name,status='unknown')
5646 read(irest2,*) totT,EK,potE,totE,t_bath
5649 ! AL 4/17/17: Now reading d_t(0,:) too
5651 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
5654 ! AL 4/17/17: Now reading d_c(0,:) too
5656 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
5659 read (irest2,*) iset
5663 end subroutine readrst
5664 !-----------------------------------------------------------------------------
5665 subroutine read_fragments
5669 use control_data, only:out1file
5672 ! implicit real*8 (a-h,o-z)
5673 ! include 'DIMENSIONS'
5677 ! include 'COMMON.SETUP'
5678 ! include 'COMMON.CHAIN'
5679 ! include 'COMMON.IOUNITS'
5680 ! include 'COMMON.MD'
5681 ! include 'COMMON.CONTROL'
5684 ! real(kind=8) :: var,ene
5686 read(inp,*) nset,nfrag,npair,nfrag_back
5688 !el from module energy
5689 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
5690 if(.not.allocated(wfrag_back)) then
5691 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5692 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5694 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
5695 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
5697 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
5698 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
5701 if(me.eq.king.or..not.out1file) &
5702 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
5703 " nfrag_back",nfrag_back
5705 read(inp,*) mset(iset)
5707 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
5709 if(me.eq.king.or..not.out1file) &
5710 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
5711 ifrag(2,i,iset), qinfrag(i,iset)
5714 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
5716 if(me.eq.king.or..not.out1file) &
5717 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
5718 ipair(2,i,iset), qinpair(i,iset)
5721 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5722 wfrag_back(3,i,iset),&
5723 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5724 if(me.eq.king.or..not.out1file) &
5725 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5726 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5730 end subroutine read_fragments
5731 !-----------------------------------------------------------------------------
5733 !-----------------------------------------------------------------------------
5737 ! implicit real*8 (a-h,o-z)
5738 ! include 'DIMENSIONS'
5739 ! include 'COMMON.CSA'
5740 ! include 'COMMON.BANK'
5741 ! include 'COMMON.IOUNITS'
5743 ! integer :: ntf,ik,iw_pdb
5744 ! real(kind=8) :: var,ene
5746 open(icsa_in,file=csa_in,status="old",err=100)
5747 read(icsa_in,*) nconf
5748 read(icsa_in,*) jstart,jend
5749 read(icsa_in,*) nstmax
5750 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5751 read(icsa_in,*) nran0,nran1,irr
5752 read(icsa_in,*) nseed
5753 read(icsa_in,*) ntotal,cut1,cut2
5754 read(icsa_in,*) estop
5755 read(icsa_in,*) icmax,irestart
5756 read(icsa_in,*) ntbankm,dele,difcut
5757 read(icsa_in,*) iref,rmscut,pnccut
5758 read(icsa_in,*) ndiff
5765 end subroutine csa_read
5766 !-----------------------------------------------------------------------------
5767 subroutine initial_write
5770 ! implicit real*8 (a-h,o-z)
5771 ! include 'DIMENSIONS'
5772 ! include 'COMMON.CSA'
5773 ! include 'COMMON.BANK'
5774 ! include 'COMMON.IOUNITS'
5776 ! integer :: ntf,ik,iw_pdb
5777 ! real(kind=8) :: var,ene
5779 open(icsa_seed,file=csa_seed,status="unknown")
5780 write(icsa_seed,*) "seed"
5782 #if defined(AIX) || defined(PGI)
5783 open(icsa_history,file=csa_history,status="unknown",&
5786 open(icsa_history,file=csa_history,status="unknown",&
5789 write(icsa_history,*) nconf
5790 write(icsa_history,*) jstart,jend
5791 write(icsa_history,*) nstmax
5792 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5793 write(icsa_history,*) nran0,nran1,irr
5794 write(icsa_history,*) nseed
5795 write(icsa_history,*) ntotal,cut1,cut2
5796 write(icsa_history,*) estop
5797 write(icsa_history,*) icmax,irestart
5798 write(icsa_history,*) ntbankm,dele,difcut
5799 write(icsa_history,*) iref,rmscut,pnccut
5800 write(icsa_history,*) ndiff
5802 write(icsa_history,*)
5805 open(icsa_bank1,file=csa_bank1,status="unknown")
5806 write(icsa_bank1,*) 0
5810 end subroutine initial_write
5811 !-----------------------------------------------------------------------------
5812 subroutine restart_write
5815 ! implicit real*8 (a-h,o-z)
5816 ! include 'DIMENSIONS'
5817 ! include 'COMMON.IOUNITS'
5818 ! include 'COMMON.CSA'
5819 ! include 'COMMON.BANK'
5821 ! integer :: ntf,ik,iw_pdb
5822 ! real(kind=8) :: var,ene
5824 #if defined(AIX) || defined(PGI)
5825 open(icsa_history,file=csa_history,position="append")
5827 open(icsa_history,file=csa_history,access="append")
5829 write(icsa_history,*)
5830 write(icsa_history,*) "This is restart"
5831 write(icsa_history,*)
5832 write(icsa_history,*) nconf
5833 write(icsa_history,*) jstart,jend
5834 write(icsa_history,*) nstmax
5835 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5836 write(icsa_history,*) nran0,nran1,irr
5837 write(icsa_history,*) nseed
5838 write(icsa_history,*) ntotal,cut1,cut2
5839 write(icsa_history,*) estop
5840 write(icsa_history,*) icmax,irestart
5841 write(icsa_history,*) ntbankm,dele,difcut
5842 write(icsa_history,*) iref,rmscut,pnccut
5843 write(icsa_history,*) ndiff
5844 write(icsa_history,*)
5845 write(icsa_history,*) "irestart is: ", irestart
5847 write(icsa_history,*)
5851 end subroutine restart_write
5852 !-----------------------------------------------------------------------------
5854 !-----------------------------------------------------------------------------
5855 subroutine write_pdb(npdb,titelloc,ee)
5857 ! implicit real*8 (a-h,o-z)
5858 ! include 'DIMENSIONS'
5859 ! include 'COMMON.IOUNITS'
5860 character(len=50) :: titelloc1
5861 character*(*) :: titelloc
5862 character(len=3) :: zahl
5863 character(len=5) :: liczba5
5865 integer :: npdb !,ilen
5869 ! real(kind=8) :: var,ene
5873 if (npdb.lt.1000) then
5874 call numstr(npdb,zahl)
5875 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5877 if (npdb.lt.10000) then
5878 write(liczba5,'(i1,i4)') 0,npdb
5880 write(liczba5,'(i5)') npdb
5882 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5884 call pdbout(ee,titelloc1,ipdb)
5887 end subroutine write_pdb
5888 !-----------------------------------------------------------------------------
5890 !-----------------------------------------------------------------------------
5891 subroutine write_thread_summary
5892 ! Thread the sequence through a database of known structures
5893 use control_data, only: refstr
5895 use energy_data, only: n_ene_comp
5897 ! implicit real*8 (a-h,o-z)
5898 ! include 'DIMENSIONS'
5900 use MPI_data !include 'COMMON.INFO'
5903 ! include 'COMMON.CONTROL'
5904 ! include 'COMMON.CHAIN'
5905 ! include 'COMMON.DBASE'
5906 ! include 'COMMON.INTERACT'
5907 ! include 'COMMON.VAR'
5908 ! include 'COMMON.THREAD'
5909 ! include 'COMMON.FFIELD'
5910 ! include 'COMMON.SBRIDGE'
5911 ! include 'COMMON.HEADER'
5912 ! include 'COMMON.NAMES'
5913 ! include 'COMMON.IOUNITS'
5914 ! include 'COMMON.TIME1'
5916 integer,dimension(maxthread) :: ip
5917 real(kind=8),dimension(0:n_ene) :: energia
5919 integer :: i,j,ii,jj,ipj,ik,kk,ist
5920 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5922 write (iout,'(30x,a/)') &
5923 ' *********** Summary threading statistics ************'
5924 write (iout,'(a)') 'Initial energies:'
5925 write (iout,'(a4,2x,a12,14a14,3a8)') &
5926 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5927 'RMSnat','NatCONT','NNCONT','RMS'
5928 ! Energy sort patterns
5933 enet=ener(n_ene-1,ip(i))
5936 if (ener(n_ene-1,ip(j)).lt.enet) then
5938 enet=ener(n_ene-1,ip(j))
5950 ist=nres_base(2,ii)+ipatt(2,i)
5952 energia(i)=ener0(kk,i)
5954 etot=ener0(n_ene_comp+1,i)
5955 rmsnat=ener0(n_ene_comp+2,i)
5956 rms=ener0(n_ene_comp+3,i)
5957 frac=ener0(n_ene_comp+4,i)
5958 frac_nn=ener0(n_ene_comp+5,i)
5961 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5962 i,str_nam(ii),ist+1,&
5963 (energia(print_order(kk)),kk=1,nprint_ene),&
5964 etot,rmsnat,frac,frac_nn,rms
5966 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5967 i,str_nam(ii),ist+1,&
5968 (energia(print_order(kk)),kk=1,nprint_ene),etot
5971 write (iout,'(//a)') 'Final energies:'
5972 write (iout,'(a4,2x,a12,17a14,3a8)') &
5973 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5974 'RMSnat','NatCONT','NNCONT','RMS'
5978 ist=nres_base(2,ii)+ipatt(2,i)
5980 energia(kk)=ener(kk,ik)
5982 etot=ener(n_ene_comp+1,i)
5983 rmsnat=ener(n_ene_comp+2,i)
5984 rms=ener(n_ene_comp+3,i)
5985 frac=ener(n_ene_comp+4,i)
5986 frac_nn=ener(n_ene_comp+5,i)
5987 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5988 i,str_nam(ii),ist+1,&
5989 (energia(print_order(kk)),kk=1,nprint_ene),&
5990 etot,rmsnat,frac,frac_nn,rms
5992 write (iout,'(/a/)') 'IEXAM array:'
5993 write (iout,'(i5)') nexcl
5995 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5997 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5998 'Max. time for threading step ',max_time_for_thread,&
5999 'Average time for threading step: ',ave_time_for_thread
6001 end subroutine write_thread_summary
6002 !-----------------------------------------------------------------------------
6003 subroutine write_stat_thread(ithread,ipattern,ist)
6005 use energy_data, only: n_ene_comp
6007 ! implicit real*8 (a-h,o-z)
6008 ! include "DIMENSIONS"
6009 ! include "COMMON.CONTROL"
6010 ! include "COMMON.IOUNITS"
6011 ! include "COMMON.THREAD"
6012 ! include "COMMON.FFIELD"
6013 ! include "COMMON.DBASE"
6014 ! include "COMMON.NAMES"
6015 real(kind=8),dimension(0:n_ene) :: energia
6017 integer :: ithread,ipattern,ist,i
6018 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
6020 #if defined(AIX) || defined(PGI)
6021 open(istat,file=statname,position='append')
6023 open(istat,file=statname,access='append')
6026 energia(i)=ener(i,ithread)
6028 etot=ener(n_ene_comp+1,ithread)
6029 rmsnat=ener(n_ene_comp+2,ithread)
6030 rms=ener(n_ene_comp+3,ithread)
6031 frac=ener(n_ene_comp+4,ithread)
6032 frac_nn=ener(n_ene_comp+5,ithread)
6033 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6034 ithread,str_nam(ipattern),ist+1,&
6035 (energia(print_order(i)),i=1,nprint_ene),&
6036 etot,rmsnat,frac,frac_nn,rms
6039 end subroutine write_stat_thread
6040 !-----------------------------------------------------------------------------
6041 subroutine readpdb_template(k)
6042 ! Read the PDB file for read_constr_homology with read2sigma
6043 ! and convert the peptide geometry into virtual-chain geometry.
6045 ! include 'DIMENSIONS'
6046 ! include 'COMMON.LOCAL'
6047 ! include 'COMMON.VAR'
6048 ! include 'COMMON.CHAIN'
6049 ! include 'COMMON.INTERACT'
6050 ! include 'COMMON.IOUNITS'
6051 ! include 'COMMON.GEO'
6052 ! include 'COMMON.NAMES'
6053 ! include 'COMMON.CONTROL'
6054 ! include 'COMMON.FRAG'
6055 ! include 'COMMON.SETUP'
6056 use compare_data, only:nhfrag,nbfrag
6057 integer :: i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity, &
6059 logical lprn /.false./,fail
6060 real(kind=8), dimension (3):: e1,e2,e3
6061 real(kind=8) :: dcj,efree_temp
6065 real(kind=8), dimension (3,20) :: sccor
6067 integer, dimension (:), allocatable :: iterter
6068 if(.not.allocated(iterter))allocate(iterter(nres))
6075 write (2,*) "UNRES_PDB",unres_pdb
6083 read (ipdbin,'(a80)',end=10) card
6084 if (card(:3).eq.'END') then
6086 else if (card(:3).eq.'TER') then
6089 itype(ires_old-1,1)=ntyp1
6090 iterter(ires_old-1)=1
6091 itype(ires_old,1)=ntyp1
6094 ! write (iout,*) "Chain ended",ires,ishift,ires_old
6097 dc(j,ires)=sccor(j,iii)
6100 call sccenter(ires,iii,sccor)
6103 ! Fish out the ATOM cards.
6104 if (index(card(1:4),'ATOM').gt.0) then
6105 read (card(12:16),*) atom
6106 ! write (iout,*) "! ",atom," !",ires
6107 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
6108 read (card(23:26),*) ires
6109 read (card(18:20),'(a3)') res
6110 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
6111 ! & " ires_old",ires_old
6112 ! write (iout,*) "ishift",ishift," ishift1",ishift1
6113 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
6114 if (ires-ishift+ishift1.ne.ires_old) then
6115 ! Calculate the CM of the preceding residue.
6119 dc(j,ires_old)=sccor(j,iii)
6122 call sccenter(ires_old,iii,sccor)
6126 ! Start new residue.
6127 if (res.eq.'Cl-' .or. res.eq.'Na+') then
6130 else if (ibeg.eq.1) then
6131 ! write (iout,*) "BEG ires",ires
6133 if (res.ne.'GLY' .and. res.ne. 'ACE') then
6137 ires=ires-ishift+ishift1
6139 ! write (iout,*) "ishift",ishift," ires",ires,
6140 ! & " ires_old",ires_old
6141 ! write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
6143 else if (ibeg.eq.2) then
6145 ishift=-ires_old+ires-1
6147 ! write (iout,*) "New chain started",ires,ishift
6150 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
6151 ires=ires-ishift+ishift1
6154 if (res.eq.'ACE' .or. res.eq.'NHE') then
6157 itype(ires,1)=rescode(ires,res,0,1)
6160 ires=ires-ishift+ishift1
6162 ! write (iout,*) "ires_old",ires_old," ires",ires
6163 ! if (card(27:27).eq."A" .or. card(27:27).eq."B") then
6166 ! write (2,*) "ires",ires," res ",res," ity",ity
6167 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
6168 res.eq.'NHE'.and.atom(:2).eq.'HN') then
6169 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
6170 ! write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
6172 write (iout,'(2i3,2x,a,3f8.3)') &
6173 ires,itype(ires,1),res,(c(j,ires),j=1,3)
6177 sccor(j,iii)=c(j,ires)
6179 if (ishift.ne.0) then
6180 ires_ca=ires+ishift-ishift1
6184 ! write (*,*) card(23:27),ires,itype(ires)
6185 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.&
6186 atom.ne.'N' .and. atom.ne.'C' .and.&
6187 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.&
6188 atom.ne.'OXT' .and. atom(:2).ne.'3H') then
6189 ! write (iout,*) "sidechain ",atom
6191 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
6195 10 if(me.eq.king.or..not.out1file) &
6196 write (iout,'(a,i5)') ' Nres: ',ires
6197 ! Calculate dummy residue coordinates inside the "chain" of a multichain
6201 ! write (iout,*) i,itype(i),itype(i+1)
6202 if (itype(i,1).eq.ntyp1.and.iterter(i).eq.1) then
6203 if (itype(i+1,1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
6204 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
6205 ! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
6206 ! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
6208 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
6209 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
6216 c(j,i)=c(j,i-1)-1.9d0*e2(j)
6220 dcj=(c(j,i-2)-c(j,i-3))/2.0
6221 if (dcj.eq.0) dcj=1.23591524223
6226 else !itype(i+1).eq.ntyp1
6228 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
6229 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
6236 c(j,i)=c(j,i+1)-1.9d0*e2(j)
6240 dcj=(c(j,i+3)-c(j,i+2))/2.0
6241 if (dcj.eq.0) dcj=1.23591524223
6246 endif !itype(i+1).eq.ntyp1
6247 endif !itype.eq.ntyp1
6249 ! Calculate the CM of the last side chain.
6252 dc(j,ires)=sccor(j,iii)
6255 call sccenter(ires,iii,sccor)
6259 if (itype(nres,1).ne.10) then
6263 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
6264 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
6271 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
6275 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
6276 if (dcj.eq.0) dcj=1.23591524223
6277 c(j,nres)=c(j,nres-1)+dcj
6278 c(j,2*nres)=c(j,nres)
6289 c(j,2*nres)=c(j,nres)
6291 if (itype(1,1).eq.ntyp1) then
6295 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
6296 call refsys(2,3,4,e1,e2,e3,fail)
6303 c(j,1)=c(j,2)-1.9d0*e2(j)
6307 dcj=(c(j,4)-c(j,3))/2.0
6313 ! Copy the coordinates to reference coordinates
6319 ! Calculate internal coordinates.
6320 if (out_template_coord) then
6321 write (iout,'(/a)') &
6322 "Cartesian coordinates of the reference structure"
6323 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
6324 "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
6326 write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)')&
6327 restyp(itype(ires,1),1),ires,(c(j,ires),j=1,3),&
6328 (c(j,ires+nres),j=1,3)
6331 ! Calculate internal coordinates.
6332 call int_from_cart(.true.,out_template_coord)
6333 call sc_loc_geom(.false.)
6335 thetaref(i)=theta(i)
6340 dc(j,i)=c(j,i+1)-c(j,i)
6341 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
6346 dc(j,i+nres)=c(j,i+nres)-c(j,i)
6347 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
6349 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
6350 ! & vbld_inv(i+nres)
6355 cref(j,i+nres,1)=c(j,i+nres)
6365 end subroutine readpdb_template
6366 !-----------------------------------------------------------------------------
6368 !-----------------------------------------------------------------------------
6369 end module io_config